NeXT Education Software Sampler 1992 Fall/NeXT Education Software Sampler 1992 Fall.iso/Mathematics/Notebooks/SigProc2.0/Packages/SignalProcessing/Support/SigProc.m

Text File | 1992-08-18 | 82.5 KB | 2,452 lines

(*  :Title:    Signal Processing Primitive Objects  *)

(*  :Authors:    Brian Evans, James McClellan  *)

(*
    :Summary:   Provide functions which are common in DSP but
        not standard in Mathematica.  Functions like
        Step and Impulse, as well as operators like
        Upsample and Shift.
 *)

(*  :Context:    SignalProcessing`Support`SigProc`  *)

(*  :PackageVersion:  2.7    *)

(*
    :Copyright:    Copyright 1989-1991 by Brian L. Evans
        Georgia Tech Research Corporation

    Permission to use, copy, modify, and distribute this software
    and its documentation for any purpose and without fee is
    hereby granted, provided that the above copyright notice
    appear in all copies and that both that copyright notice and
    this permission notice appear in supporting documentation,
    and that the name of the Georgia Tech Research Corporation,
    Georgia Tech, or Georgia Institute of Technology not be used
    in advertising or publicity pertaining to distribution of the
    software without specific, written prior permission.  Georgia
    Tech makes no representations about the suitability of this
    software for any purpose.  It is provided "as is" without
    express or implied warranty.
 *)

(*  :History:    signal processing, operators, functions, symbolic processing *)

(*  :Keywords:    *)

(*
    :Source:         Oppenheim, A., and Schafer, R.  {Discrete-Time
        Signal Processing}.  1989.

             Myers, C.  {Signal Representation for Symbolic and
        Numeric Processing}.  MIT Ph.D. Thesis.  1986.

             Covell, M.  {An Algorithm Design Environment}.
        MIT Ph. D. Thesis, 1989.

             Bamberger, R. {The Directional Filter Bank:  A Multirate
        Filter Bank for the Directional Decompostion of Images},
        Georgia Tech Ph. D. Thesis, 1990.
 *)

(*  :Warning:    *)

(*  :Mathematica Version:  1.2 or 2.0  *)

(*  :Limitation:  *)

(*
    :Discussion:    This package introduces the functions (signals) and
          operators (systems) common in signal processing but missing
          in Mathematica.
            Thanks to John Mitchell at the University of Chicago
          for writing patches for the Summation operator and the Delta
          function.
*)

(*  :Functions: Aliasby            operator
        AliasSinc        alias
        ASinc            alias
        CConvolve        operator
        CircularShift        operator
        Convolve        operator
        Continuous        option value
        CPulse            function
        CStep            function
        Delta            function
        DeltaPlot        function
        DFT            operator
        DTFT            operator
        Discrete        option value
        Difference        operator
        Dirichlet        function
        DiscreteGraphics    routine
        Domain            option
        Downsample        operator
        DummyVariables        function
        Extent1D        function
        FT            operator
        GetDeltaFunctions    function
        Impulse            function
        Interleave        operator
        IntervalsToFunction    routine
        InvDFT            operator
        InvFT            operator
        InvL            operator
        InvZ            operator
        L            operator
        LineImpulse        function
        OperatorInOperExpr    routine
        OperatorName        routine
        ParametersInOperExpr    routine
        PolyphaseDownsample    operator
        PolyphaseUpsample    operator
        Periodic        operator
        Pulse            function
        RationalGCD        routine
        Rev            operator
        RewriteSummations    routine
        ScaleAxis        operator
        ScalingFactor        routine
        SequenceToFunction    routine
        SequencePlot        routine
        Shift            operator
        SignalCleanup        routine
        SignalPlot        function
        SignalsInOperExpr    routine
        Sinc            function
        SPSimplify        function
        Step            function
        Summation        operator
        TheFunction        routine
        ToContinuous        routine
        ToDiscrete        routine
        Unit            function
        Upsample        operator
        UpsampledFunction    routine
        UpsampleFactor        routine
        UpsampleSequence    routine
        Z            operator

    This package provides signal primitives (functions), system
    primitives (operators), and supporting routines.   Each new
    object is placed into the SPsignals, SPsystems, or SPfunctions
    list accordingly (at end of file).

    Each signal primitive (function) is specified as follows:

    (1) default arguments (if any),
    (2) definition (so that it behaves as a function),
    (3) multidimensional separability (if any),
    (4) automatic simplification rules,
    (5) manual simplification rules (augment Simplify),
    (6) output form rules (including TeX form rules), and
    (7) derivative and integration rules (if any)

    Each system primitive (operator) is specified as follows:

    (1) multidimensional separability (if any),
    (2) manual simplification rules (augment Simplify),
    (3) definition (so that it can be converted into a function form),
    (4) output form rules (including TeX form rules), and
    (5) specify which parameter is (are) the variable(s)
 *)



(*  B E G I N     P A C K A G E  *)

BeginPackage[ "SignalProcessing`Support`SigProc`",
          "SignalProcessing`Support`SupCode`" ]


If [ TrueQ[ $VersionNumber >= 2.0 ],
     Off[ General::spell ];
     Off[ General::spell1 ] ]


(*  U S A G E     I N F O R M A T I O N  *)

$DeltaFunctionScaling::usage =
    "The variable $DeltaFunctionScaling dictates how to scale the heighth \
    of Delta functions for the purposes of plotting:  None or Scaled."

Aliasby::usage =
    "Aliasby[m,w][x] represents the aliasing of x, a continuous-frequency \
    signal, by replicating it every (2 Pi / m) in w and dividing the \
    result by m. \
    It is always rewritten in terms of the Periodic operator."

CConvolve::usage =
    "CConvolve[t][x1, x2, ...] represents the continuous convolution in t \
    of functions x1, x2, .... \
    CConvolve[{t1, t2, ...}][x1, x2, ...] represents the multidimensional \
    continuous convolution in t1, t2, ... of functions x1, x2, ..."

CircularShift::usage =
    "CircularShift[n0, N, n] represents the circular shifting of a \
    sequence of length N in index variable n by n0. \
    That is, the sequence is treated as a ring buffer of length N \
    so each sample at index n is shifted to the right by
    (n + n0) mod N places. \
    See also Shift and RotateRight."

Convolve::usage =
    "Convolve[n][x1, x2, ...] represents the discrete convolution in n of \
    functions x1, x2, .... \
    Convolve[{n1, n2, ...}][x1, x2, ...] represents the multidimensional \
    discrete convolution in n1, n2, ... of functions x1, x2, ..."

Continuous::usage =
    "Continuous is a possible value for a signal's domain."

CPulse::usage =
    "CPulse[l,t] defines a pulse which begins at t=0 and ends at t = l. \
    The CPulse has value 1 within the range (0,l), \
    0 outside this range, and 1/2 at the points t=0 and t=l. \
    A continuous-pulse center at the origin is written as \
    CPulse[l, t + l/2] or Shift[-l/2, t][CPulse[l, t]]."

CStep::usage =
    "CStep[t], a.k.a. Unit[-1][t], is the unit step function \
    which is 1 for t > 0, 0 for t < 0, and 1/2 at t = 0. \
    It is commonly used for continuous expressions t. \
    See also Step and Unit."

Delta::usage =
    "Delta[expr] is the Dirac delta function. \
    The area under this functions is 1 but it only has value \
    at the origin. \
    That is, Integrate[ Delta[t] g[t], {t, t1, t2} ] is g[0] \
    if t1 <= 0 <= t2, 0 otherwise. \
    It differs from the Kronecker delta function Impulse[t]."

DeltaIntegrate::usage =
    "DeltaIntegrate[fun, {x, xlower, xupper}] and DeltaIntegrate[fun, x] \
    are special hooks for the built-in Integrate routine."

DeltaPlot::usage =
    "DeltaPlot is a object that can generate a plot of only Dirac delta \
    functions or a combination of Dirac delta functions and another plot. \
    In the first case, either four or six arguments are passed: \
    deltafuns, t, tmin, tmax, ymin, ymax.\
    In the second case, either five or seven arguments are passed: \
    deltafuns, t, tmin, tmax, plot, ymin, ymax. \
    In both cases, ymin and ymax are the optional arguments. \
    DeltaPlot returns a list of plots suitable for sending to Show. \
    Note that 1-D Delta functions are plotted as arrows. \
    The variable $DeltaFunctionScaling dictates how Delta functions \
    are scaled when plotted:  None or Scaled."

DFT::usage =
    "DFT[L, n, k] is the forward discrete Fourier transform operator. \
    Here, L is the DFT length(s), n is the time-domain variable(s), and \
    k is the frequency-domain variable(s). \
    Multidimensional DFT's are supported. \
    Applying TheFunction to the DFT operator will invoke \
    the DFT rule base (object DFTransform) if loaded."

DTFT::usage =
    "DTFT[n, w] is the forward discrete-time Fourier transform operator. \
    Here, n is the time-domain variable(s), and w is the frequency-domain \
    variable(s). \
    Multidimensional DTFT's are supported. \
    Applying TheFunction to the DTFT operator will invoke the DTFT \
    rule base (object DTFTransform) if loaded."

Difference::usage =
    "Difference[k,n] represents the k-th backward difference with respect \
    to the variable n. \
    The first-order difference of f equals f[n] - f[n-1], \
    the second-order difference is f[n] - 2 f[n-1] + f[n-2], and so on. \
    Applying this operator to the sequence x[n] would be written as \
    Difference[k,n][ x[n] ]."

Dirichlet::usage =
    "Dirichlet[N, w] is the aliased sinc function. \
    which is defined as Sin[N w / 2] / ( N Sin [w / 2] ). \
    Aliases for Dirichlet are ASinc, AliasSinc, and AliasedSinc."

Discrete::usage =
    "Discrete is a possible value for a signal's domain."

DiscreteGraphics::usage =
    "DiscreteGraphics[{{x1,y1}, {x2,y2}, ... {xN,yN}}] will plot \
    the points as lines drawn from the x-axis to the y-value. \
    It takes two optional arguments:  a line thickness and \
    a point thickness. \
    Each is a real value from 0 to 1 indicating a width that \
    is a fraction of the size of the entire graph."

Domain::usage =
    "Domain is an option for many signal processing functions. \
    It is either Discrete or Continuous."

Downsample::usage =
    "Downsample[m,n] represents the downsampling operator. \
    Downsampling resamples a function at the points n = 0, m, 2m, 3m, ... \
    In N dimensions, n is a list of N indices and m is an N x N matrix. \
    Applying TheFunction to f[n] which is downsampled in one dimension, \
    (written as Downsample[m, n][ f[n] ]) will produce f[m n]."

DummyVariables::usage =
    "DummyVariables[dimension, variable] generates dimension dummy \
    variables with variable as the root. \
    For example, n can be generated by DummyVariables[1, n] or \
    DummyVariables[0, n]. \
    As another example, DummyVariables[3, t] returns {t1, t2, t3}. \
    In this cases, the generated variables have the same context \
    as variable. \
    DummyVariables[dimension, variable, context] can \
    be used to specify another context/package."

Extent1D::usage =
    "Extent1D[expression, var] returns an interval \
    {lower-limit, upper-limit} outside of which the function is zero. \
    This routine assumes that var is a real variable."

FT::usage =
    "FT[t, w] is the forward continuous Fourier transform operator. \
    Here, t is the time-domain variable(s) and w is the \
    frequency-domain variable(s). \
    Applying TheFunction to this operator will invoke the continuous-time \
    Fourier transform rule base (object CTFTransform) if loaded."

GetDeltaFunctions::usage =
    "GetDeltaFunctions[function, variable(s)] will extract \
    all delta functions from the argument function.  See also DeltaPlot."

Impulse::usage =
    "Impulse[n] is the unit Kronecker Delta function. \
    At n = 0, the function evaluates to 1. \
    Elsewhere, it evaluates to 0."

Interleave::usage =
    "Interleave[n][x0, x1, ...] interleaves samples of signals \
    x0, x1, ..., which are functions of n. \
    This is only applicable to discrete signals."

IntervalsToFunction::usage =
    "IntervalsToFunction[interlist, n, step, pulse] \
    will translate the list of functions and intervals (interlist) \
    into a signal expression as a function of n. \
    Each element in the interlist is a list of three elements: \
    <function_of_n>, <n_minus>, and <n_plus>. \
    The function <function_of_n> will be defined for \
    <n_minus> <= n <= <n_plus>. \
    Note that <n_minus> can be an integer or -Infinity and \
    that <n_plus> can be an integer or Infinity. \
    The default values are n = Global`n, \
    step = Step, and pulse = Pulse. \
    This function also works for when n represents a continuous variable \
    (step = CStep and pulse = CPulse)."

InvDFT::usage =
    "InvDFT[L, k, n] is the inverse discrete Fourier transform operator. \
    Here, L is the DFT length(s), n is the time-domain variable(s), and \
    k is the frequency-domain variable(s). \
    Multidimensional inverse DFT's are supported. \
    Applying TheFunction to the InvDFT operator will invoke the \
    inverse DFT rule base (object InvDFTransform) if loaded."

InvDTFT::usage =
    "InvDTFT[w, n] is the inverse discrete-time Fourier transform \
    operator. \
    Here, n is the time-domain variable(s), and w is the \
    frequency-domain variable(s). \
    Multidimensional inverse DTFT's are supported. \
    Applying TheFunction to the InvDTFT operator will invoke the \
    inverse DTFT rule base InvDTFTransform if loaded."

InvFT::usage =
    "InvFT[w, t] is the inverse discrete Fourier transform operator. \
    Here, t is the time-domain variable(s) and w is the \
    frequency-domain variable(s). \
    Multidimensional Fourier transforms are supported. \
    Applying TheFunction to this operator will invoke the inverse \
    continuous-time Fourier transform rule base InvCTFTransform if loaded."

InvL::usage =
    "InvL[s, t] is the inverse Laplace transform operator. \
    Applying TheFunction to the InvL operator will invoke the \
    bilateral inverse Laplace rule base InvLaPlace if loaded."

InvZ::usage =
    "InvZ[z, n] is the inverse z-transform operator. \
    Applying TheFunction to the InvZ operator will invoke the inverse \
    z-transform rule base InvZTransform if loaded."

L::usage =
    "L[t, s] is the forward Laplace transform operator. \
    Applying TheFunction to the L operator will invoke the bilateral \
    Laplace rule base if loaded."

LineImpulse::usage =
    "LineImpulse[nlist, coefflist] represents a line impulse where \
    nlist is a list of variables like {n1,n2,n3} and coefflist is \
    a corresponding list of coefficients like {1,2,3}. \
    For the previous lists, the line impulse is a set of impulse \
    along the line n1 = 2 n2 = 3 n3."

OperatorInOperExpr::usage =
    "OperatorInOperExpr[ operator_expression ] returns the full \
    operator in operator_expression. \
    For example, Shift[2,n] would be returned from \
    OperatorInOperExpr[ Shift[2,n][x[n]] ]."

OperatorName::usage =
    "OperatorName[ operator_expression ] returns the head of the \
    expression. \
    For example, OperatorName[ Shift[2,n][x[n]] ] will return Shift."

ParametersInOperExpr::usage =
    "ParametersInOperExpr[ operator_expression ] returns the parameters \
    of the operator in operator_expression. \
    For example, (2,n) would be returned from \
    ParametersInOperExpr[ Shift[2,n][x[n]] ]."

Periodic::usage =
    "Periodic[k,n][f] indicates that f is a function of n which is \
    periodic with period of k."

PolyphaseDownsample::usage =
    "PolyphaseDownsample[m, h, n][x] is equivalent to downsampling \
    x by m (w/r to n) and then filtering (w/r to n) by a polyphase \
    form of the FIR filter h. \
    Hence, the filtering is carried out at the lower sampling rate."

PolyphaseUpsample::usage =
    "PolyphaseUpsample[l, h, n][x] is equivalent to filtering (w/r to n) \
    x by a polyphase form of the FIR filter h and then upsampling the \
    result by l (w/r to n). \
    Hence, the filtering is carried out at the lower sampling rate."

Pulse::usage =
    "Pulse[len] and Pulse[len, n] define a pulse which begins at \
    n = 0 and ends at n = len - 1. \
    A discrete pulse centered at the origin is written as \
    Pulse[len, n + (len - 1)/2] or Shift[-(len - 1)/2, n][Pulse[len, n]]. \
    Sometime, a pulse will be automatically simplified; \
    e.g., Pulse[0, n] is zero and Pulse[1, n] is really Impulse[n]."

RationalGCD::usage =
    "RationalGCD[e1, e2, ...] returns the greatest common divisor of the \
    numerators of elements ei divided by the greatest common divisor of \
    the denominators of elements ei. \
    Every element must be a polynomial or a number. \
    Each real-valued and complex-valued element is rationalized \
    before GCD's are computed."

Rev::usage =
    "Rev[n][x] represents the operation of reversing the direction of x, \
    with respect to the variables(s) n."

RewriteSummations::usage =
    "RewriteSummations[f, t, tlower, tupper] rewrites Summation \
    and Periodic operators that appear in f so that the domain \
    of the new expression includes t in [tlower, tupper]. \
    It relies on the Extent1D function to compute the domain \
    of the arguments to the Periodic and Summation operators."

ScaleAxis::usage =
    "ScaleAxis[l,w][x] represents the compression by a factor of l (an \
    integer) of the continuous w axis of x, which is usually a \
    continuous-frequency signal."

ScalingFactor::usage =
    "ScalingFactor[expr, variable] returns the greatest common divisor \
    of all coefficients of powers of z in expr. \
    The powers of z that are functions of z are not considered. \
    For example, ScalingFactor[ a z + a^2 z^2, z ] returns a. \
    This function enables the implementation of the similarity \
    property for transforms."

SequencePlot::usage =
    "SequencePlot[f, {n, start, end}] plots real-valued 1-D sequences \
    and has the same options as Show. \
    You will have to apply Re or Im to f to plot the real or \
    imaginary values of the sequence, respectively. \
    SequencePlot[f, {n1, start1, end1}, {n2, start2, end2}] plots \
    2-D sequences and supports the same options as Plot3D."

SequenceToFunction::usage =
    "SequenceToFunction[seq], SequenceToFunction[seq, n], and \
    SequenceToFunction[seq, n, noffset] return the sequence seq, \
    {x[noffset], ..., x[noffset + N - 1]}, in symbolic form as a sum of \
    impulse functions. \
    For example, SequenceToFunction[{1, 2, 3}] -> \
    Impulse[n] + 2 Impulse[n - 1] + 3 Impulse[n - 2]. \
    For multidimensional sequences, the argument n must be a list of \
    variables of appropriate length. \
    For example, SequenceToFunction[{{1,2},{3,4}}, {n1, n2}] -> \
    Impulse[n1,n2] + 2 Impulse[n1,n2-1] + 3 Impulse[n1-1,n2] + \
    4 Impulse[n1-1, n2-2]."

Shift::usage =
    "Shift[m,v][x] represents the shifting of the expression x by m \
    samples to the right in the v direction, i.e. x[v] --> x[v - m]."

SignalCleanup::usage =
    "SignalCleanup[ expr, t, nstart, nend ] converts the signal \
    expression expr in variables t into a normal Mathematica expression \
    that can be plotted over the interval [nstart, nend]. \
    It returns a list of two expressions:  the first consists of the
    Delta functions extracted from expr, and the second is the rest
    of the expression."

SignalPlot::usage =
    "SignalPlot[f, {t, start, end}] will plot f(t) as an one-dimensional, \
    continuous-time function. \
    It will show the real part as solid lines, \
    and the imaginary part as dashed lines. \
    Delta functions are plotted as upward pointing arrows. \
    SignalPlot[f, {t1, start1, end1}, {t2, start2, end2}] treats f \
    as a function of two variables t1 and t2. \
    SignalPlot supports the same options as Plot for 1-D signals \
    (functions) and Plot3D for 2-D signals (functions)."

SignalsInOperExpr::usage =
    "SignalsInOperExpr[ operator_expression ] returns the signal(s) \
    of the operator in operator_expression. \
    For example, { x[n] } would be returned from \
    SignalsInOperExpr[ Shift[2,n][x[n]] ]. \
    Another version SignalsInOperExpr[ operator_expression, operator ] \
    behaves as the first version if operator is the head of \
    operator_expression; otherwise, operator_expression is returned."

Sinc::usage =
    "Sinc[t] is the unit sinc function which is the ratio of \
    Sin[t] over t.  Note that Sinc[0] = 1."

SPSimplify::usage =
    "SPSimplify[e] will simplify expression e according to Mathematica's \
    simplification rules and the rules in SPSimplificationRules. \
    At each iteration, Simplify is called and SPSimplificationRules \
    is applied to the result. \
    Unlike Simplify, SPSimplify supports a Dialogue option. \
    You can specify which symbols to treat as real variables by using \
    the Variables option."

Step::usage =
    "Step[n] is the unit step function which is 1 for n >= 0 \
    and 0 otherwise. \
    It is commonly used for discrete n. \
    This function differs from CStep only at the origin."

Summation::usage =
    "Summation[i][f], Summation[i, iend][f], and \
    Summation[i, istart, iend, inc][f] represents the summation of f \
    with respect to variable i:  i = (istart) to (iend) step (inc). \
    The Summation operator is an abstraction of Mathematica's Sum object. \
    Applying the object TheFunction to a Summation operator will invoke \
    the Sum object if istart, iend, and inc are numbers."

TheFunction::usage =
    "TheFunction[object] returns the mathematical function embedded \
    in object. \
    For list and transform objects, this function returns \
    the first element after the head of the object; otherwise, the \
    contents of the slot TheFunction are returned. \
    For signal processing operations, \
    the equivalent mathematical operation is returned."

ToContinuous::usage =
    "ToContinuous[expr] replaces any discrete operator or signal \
    in expr with its continuous equivalent."

ToDiscrete::usage =
    "ToDiscrete[expr] replaces any continuous operator or signal \
    in expr with its discrete equivalent."

Unit::usage =
    "Unit[n][t] is the nth order continuous unit function in t. \
    Unit[0][t] is the Dirac delta function (see Delta). \
    Unit[-1][t] is the continuous step function (see CStep)."

Upsample::usage =
    "Upsample[k,n][f] represents the upsampling operation on f, \
    which is either an expression of n or a sequence (list) of values. \
    Upsampling inserts k-1 zeroes between every two samples. \
    In N dimensions, n is a list of N indices and k is an N x N matrix."

UpsampledFunction::usage =
    "UpsampledFunction[m, fill, n, upf] represents the upsampling of f. \
    In one-dimension, the upsampled form of f, upf, \
    is obtained by evaluating n -> n / m."

UpsampleFactor::usage =
    "UpsampleFactor[f, n] returns the upsampling index which can be \
    a negative or positive number. \
    If this routine returns 1 or -1, then no upsampling has occurred; \
    if it returns 0, then f is constant with respect to n."

UpsampleSequence::usage =
    "UpsampleSequence[seq,m,fill] inserts m-1 fill values between \
    every two elements of seq."

Z::usage =
    "Z[n, z] is the forward z-transform operator. \
    Applying TheFunction to the Z operator will invoke the z-transform \
    rule base ZTransform if loaded."

(*  E N D     U S A G E     I N F O R M A T I O N  *)


Begin["`Private`"]


(*  C O N S T A N T S  *)

axesdefaultvalue = If [ TrueQ[$VersionNumber >= 2.0], True, Automatic]

$DeltaFunctionScaling = Scaled


(*  M E S S A G E S  *)

CircularShift::badlength =
    "The sequence length parameter of the CircularShift \
    operator must be non-zero."

DeltaPlot::unresolved = "These symbols were assigned a value of 1: ``."

Downsample::convert =
    "Converting downsampling vector into a diagonal matrix: ``."
Downsample::badmD =
    "Downsampling factor has different dimensions than downsampled variables."

Pulse::notsym = "Last argument to Pulse is not a symbol."

Shift::badmD =
    "Shift vector has different dimensions than the specified variables."

SequencePlot::nolist =
    "It does not make sense to plot a list of 2-D sequences."

SignalPlot::nolist =
    "It does not make sense to plot a list of 2-D functions."

Summation::extent =
    "Could not compute which part of the summation should be plotted."

Upsample::convert =
    "Converting upsampling vector into a diagonal matrix: ``."
Upsample::badmD =
    "Upsampling factor has different dimensions than upsampled variables."



(*  L O C A L     F U N C T I O N S  *)

(*  breakup  *)
breakup[a_, z_] := a /; breakupQ[z][a]
breakup[a_Times, z_] := Select[a, breakupQ[z]]
breakup[x_, z_] := 1

breakupQ[z_][a_] := FreeQ[a, z]

(*  discBreakup  *)
discBreakup[a_, z_] := a /; discBreakupQ[z][a]
discBreakup[a_Times, z_] := Select[a, discBreakupQ[z]]
discBreakup[x_, z_] := 1

discBreakupQ[z_][a_] := FreeQ[a, z] && Implies[ NumberQ[a], IntegerQ[a] ]

(*  gcdfilter --  rationalize floating point numbers        *)
gcdfilter[x_Complex] := ToCollection[ gcdfilter[Re[x]], gcdfilter[Im[x]] ]
gcdfilter[x_Real] := Rationalize[x]
gcdfilter[x_] := x

(*  informatargs and informat  *)
SetAttributes[informatargs, {Listable}]
informatargs[m_, n_] := ToString[ StringForm["``  in  ``", m, n] ]

informat[op_, m_, n_] :=
Format[ SequenceForm[op, Subscript[ ColumnForm[ToList[informatargs[m, n]]] ]] ]

(* intdivide --  kludge for Quotient primitive which has a problem  *)
(*               under Mma 1.2 when the second argument is negative *)
If [ TrueQ[ $VersionNumber >= 2.0 ],
     intdivide[n_, m_] := Quotient[n, m],
     intdivide[n_, m_] := Sign[m] Quotient[n, m Sign[m]] ]

(*  inTeXformat  *)
inTeXformat[op_, m_, n_Symbol] := Subscripted[ op[m, n] ]
inTeXformat[op_, m_, n_List] := Subscripted[ op[texFix[m], texFix[n]] ]

(*  mymax  *)
mymax[x_, y_] := minmaxSimplify[ Max[Expand[x], Expand[y]] ]

(*  mymin  *)
mymin[x_, y_] := minmaxSimplify[ Min[Expand[x], Expand[y]] ]

(*  minmaxSimplify  *)
minmaxSimplify[x_] :=
    x //. SignalProcessing`Support`SupCode`Private`MinMaxRules

(*  texFix  *)
texFix[a_?VectorQ] := MatrixForm[ Transpose[{a}] ]
texFix[a_] := delimitedMatrix[ a ]

(*  delimitedMatrix  *)
delimitedMatrix/: Format[ delimitedMatrix[a_], TeXForm ] :=
    StringForm["\\left[ `` \\right]", MatrixForm[a] ]


(*  D E T E R M I N A T I O N     O F      V A R I A B L E S    *)

DummyVariables[i_, v_] :=
    If [ i < 2,
         v,
         Table[GenerateSymbol[v, index], {index, 1, i}] ] /;
    IntegerQ[i] && ( i >= 0 )

DummyVariables[i_, v_, context_] :=
    If [ i < 2,
         GenerateSymbol[v, " ", context],
         Table[GenerateSymbol[v, index, context], {index, 1, i}] ] /;
    IntegerQ[i] && ( i >= 0 )


(*  F U N C T I O N S  *)

(*  Aliasby  *)
Aliasby/: Aliasby[m_, w_][f_] := 1/m Periodic[2 Pi / m, w][ f ]

(*  old definitions  *)
(*
Aliasby/: Aliasby[1, w_][x_] := x

Aliasby/: Format[ Aliasby[m_, w_] ] := informat[ Aliasby, m, w ]
Aliasby/: Format[ Aliasby[m_, w_], TeXForm ] := inTeXformat[ Aliasby, m, w ]

Aliasby/: GetOperatorVariables[ Aliasby[m_, w_] ] := w
 *)

(*  CConvolve --  continuous convolution  *)
CConvolve/: CConvolve[{t_Symbol}] := CConvolve[t]
CConvolve/: CConvolve[t1_, trest__] := CConvolve[{t1, trest}]

CConvolve/: CConvolve[t_Symbol] [x_, c_. Delta[t_ + m_.]] :=
    (c x) /. t -> (-m) /; FreeQ[m, t]
CConvolve/: CConvolve[t_Symbol] [c_. Delta[t_ + m_.], x_] :=
    (c x) /. t -> (-m) /; FreeQ[m, t]

CConvolve/: Extent1D[ CConvolve[t_][x__], t_ ] :=
    Apply[ Plus, Map[ Extent1D[t], {x} ] ]

CConvolve/: CConvolve[tlist_List] [x_, c_. Delta[t_ + m_.]] :=
    CConvolve[ Complement[tlist, {t}] ][x /. t -> (-m), c] /;
    MemberQ[tlist, t] && FreeQ[m, t]
CConvolve/: CConvolve[tlist_List] [c_. Delta[t_ + m_.], x_] :=
    CConvolve[ Complement[tlist, {t}] ][c, x /. t -> (-m)] /;
    MemberQ[tlist, t] && FreeQ[m, t]

CConvolve/: CConvolve[t_][x_, z_?ZeroQ] := z
CConvolve/: CConvolve[t_][z_?ZeroQ, x_] := z

CConvolve/: Format[CConvolve[t_]] :=
        Format[ Subscripted[CConvolve[t]] ]
CConvolve/: Format[CConvolve[t_][x1_, x2_], TeXForm] :=
        StringForm[ "`` \\star_{``} ``", Format[x1], t, Format[x2] ]
CConvolve/: Format[CConvolve[t_][x1_, x2_, rest__], TeXForm] :=
        Block [    {format, i, length, operstr, xlist},
            xlist = {x1, x2, rest};
            length = Length[xlist];
            operstr = ToString[StringForm["\\star_{``}", t]];
            str = StringJoin[" ", operstr, " ``"];
            format = Apply[ StringJoin,
                    Table[ str, {i, 1, length} ] ];
            format = StringJoin["`` ", operstr, " ``", format];
            Apply[StringForm, {format, Map[Format, xlist]}] ]

(*  CircularShift  *)
CircularShift/: CircularShift[n0_, N_] := CircularShift[n0, N, Global`n]

CircularShift/: CircularShift[n0_, 0, n_][x_] :=
    MyMessage[CircularShift::badlength, x]

CircularShift/: Extent1D[ CircularShift[n0_, N_, n_][f_], n_ ] := { 0, N-1 }

CircularShift/: TheFunction[ CircularShift[n0_Integer, N_Integer, n_][f_List] ] :=
    RotateRight[ Take[f, N], n0 ]

CircularShift/: TheFunction[ CircularShift[n0_, N_, n_][f_] ] :=
    TheFunction[f Pulse[N, n]] /. ( n -> Mod[n + n0, N] )

CircularShift/: N[ CircularShift[n0_, N_, n_][f_] ] :=
    N[ TheFunction[CircularShift[n0, N, n][f]] /. ( n -> Mod[n + n0, N] ) ]

CircularShift/: Format[ CircularShift[n0_, N_, n_] ] :=
    informat[CircularShift, Mod[n + n0, N], n]
CircularShift/: Format[ CircularShift[m_,n_], TeXForm ] :=
    inTeXformat[CircularShift, Mod[n + n0, N], n]

(*  CPulse  *)
CPulse[l_] := CPulse[l, Global`t]            (* default args     *)

CPulse/: CPulse[l_, t_] := 1     /; N[0 < t < l]    (* definition        *)
CPulse/: CPulse[l_, t_] := 0     /; N[(t < 0) || (t > l)]
CPulse/: CPulse[l_, t_] := 1/2   /; N[(t == 0) || (t == l)]
CPulse/: CPulse[Infinity, t_] := CStep[t]
CPulse/: Extent1D[ CPulse[l_, a_. t_ + b_.], t_ ] :=
        Extent1DOrder[-b/a, (l-b)/a] /;
        FreeQ[{a,b,l}, t]

CPulse/: Format[ CPulse[l_, n_] ] :=            (* output forms     *)
        Format[ CPulse Subscript[l] [n] ]

CPulse/: Format[ CPulse, TeXForm ] := "\\sqcap"
CPulse/: Format[ CPulse[l_, n_], TeXForm ] :=
        StringForm[ "\\sqcap_{``} (``)", l, n ]

Derivative[0,1][CPulse] :=                (* derivative rule *)
    ( Delta[#2] + Delta[#2 - #1] )&

Unprotect[Integrate]                    (* integration rules *)
Integrate/: Integrate[f_. CPulse[L_, a_. x_ + b_.] + c_.,
        vars1___, x_Symbol, vars2___] :=
    1/a CPulse[a x + b] ( Integrate[f, vars1, x, vars2] /. x -> a x + b ) +
      Integrate[c, vars1, x, vars2] /;
    FreeQ[{a, b, L}, x]
Integrate/: Integrate[expr_. CPulse[L_, a_. x_ + b_.] + c_.,
            vars1___, {x_Symbol, x1_, x2_}, vars2___] :=
    Sign[x2 - x1] *
    cpulseintersect[ Min[x1, x2], Max[x1, x2],
             Min[-b/a, (L - b)/a], Max[-b/a, (L - b)/a] ] *
    Integrate[expr,
          vars1,
          {x, Max[Min[x1,x2], Min[-b/a, (L - b)/a]],
              Min[Max[x1,x2], Max[-b/a, (L - b)/a]]},
              vars2] +
    Integrate[c, vars1, {x, x1, x2}, vars2] /;
    FreeQ[{a, b, L}, x]
Protect[Integrate]

cpulseintersect[x1_, x2_, x3_, x4_] := Step[x4 - x1] Step[x2 - x3]

(* CStep  *)
CStep/: CStep[Infinity] := 1                (* definition         *)
CStep/: CStep[t_] := 1        /; N[ t > 0 ]
CStep/: CStep[t_] := 0        /; N[ t < 0 ]
CStep/: CStep[0] := 1/2
CStep/: CStep[0.] := 0.5
CStep/: CStep[-Infinity] := 0

CStep/: Extent1D[ CStep[a_. t_ + b_.], t_ ] :=
        Extent1DOrder[-b/a, Sign[a] Infinity] /;
        FreeQ[{a,b}, t]

CStep/: CStep[t_, args__] := CStep[t] CStep[args]    (* multidimensional *)
CStep/: CStep[List[args__]] := CStep[args]        (* separability     *)

CStep/: CStep[n_Symbol + m_.] CStep[n_Symbol  + k_.] :=    (* automatic        *)
    CStep[n + Min[m,k]] /;                (* simplification   *)
    NumberQ[m] && NumberQ[k]            (* rules        *)
CStep/: CStep[n_Symbol + m_.] CStep[-n_Symbol + k_.] :=
    CPulse[k + m, n + m] /;
    N[ -m <= k ]
CStep/: CStep[n_Symbol + m_.] CStep[-n_Symbol + k_.] := 0 /;
    N[ -m > k ]

CStep/: Simplify[CStep[a_ n_Symbol + b_.]] :=        (* man. simp. rules *)
    CStep[n + b/a] /;
    MyFreeQ[{a,b}, n] && Implies[ NumberQ[N[a]], N[a > 0] ]
CStep/: Simplify[CStep[- a_ n_Symbol + b_.]] :=
    CStep[-n + b/a] /;
    MyFreeQ[{a,b}, n] && Implies[ NumberQ[N[a]], N[a > 0] ]

CStep/: Format[ CStep, TeXForm ] := "u_{-1}"        (* output form     *)

Derivative[1][CStep] := Delta                (* derivative rule *)

Unprotect[Integrate]                    (* integration rules *)
Integrate/: Integrate[f_. CStep[a_. x_ + b_.] + c_.,
        vars1___, x_Symbol, vars2___] :=
    1/a CStep[a x + b] ( Integrate[f, vars1, x, vars2] /. x -> a x + b ) +
      Integrate[c, vars1, x, vars2] /;
    FreeQ[{a, b}, x]
Integrate/: Integrate[expr_. CStep[a_. x_ + b_.] + c_.,
            vars1___, {x_Symbol, x1_, x2_}, vars2___] :=
    Sign[x2 - x1] *
    cstepintersect[x1, x2, a, b] *
    Integrate[expr,
          vars1,
          {x, csteplbound[x1, x2, a, b], cstepubound[x1, x2, a, b]},
              vars2] +
    Integrate[c, vars1, {x, x1, x2}, vars2] /;
    FreeQ[{a, b}, x]
Protect[Integrate]

cstepintersect[x1_, x2_, a_, b_] :=
    If [ Positive[a],
         Step[Max[x1, x2] + b/a],
             Step[-b/a - Min[x1,x2]] ]

csteplbound[x1_, x2_, a_, b_] :=
    minmaxSimplify[Max[Min[x1, x2], -b/a]] /; Positive[a]
csteplbound[x1_, x2_, a_, b_] :=
    minmaxSimplify[Min[x1, x2]] /; Negative[a]
cstepubound[x1_, x2_, a_, b_] :=
    minmaxSimplify[Max[x1, x2]] /; Positive[a]
cstepubound[x1_, x2_, a_, b_] :=
    minmaxSimplify[Min[Max[x1, x2], -b/a]] /; Negative[a]

Format[csteplbound[x1_, x2_, a_, b_]] :=
    SequenceForm["LeftEndpointOf[", {x1,x2}, " and step at ", -b/a, "]"]
Format[cstepubound[x1_, x2_, a_, b_]] :=
    SequenceForm["RightEndpointOf[", {x1,x2}, " and step at ", -b/a, "]"]

(*  Convolve --  discrete convolution  *)
Convolve/: Convolve[{n_Symbol}] := Convolve[n]
Convolve/: Convolve[n1_, nrest__] := Convolve[{n1, nrest}]

Convolve/: Convolve[n_] [x_, c_. Impulse[n_ + m_.]] :=
    (c x) /. n -> (-m) /; FreeQ[m, t]
Convolve/: Convolve[n_] [c_. Impulse[n_ + m_.], x_] :=
    (c x) /. n -> (-m) /; FreeQ[m, t]

Convolve/: Extent1D[ Convolve[n_][x__], n_ ] :=
    Apply[ Plus, Map[ Extent1D[n], {x} ] ]

Convolve/: Convolve[nlist_List] [x_, c_. Impulse[n_ + m_.]] :=
    Convolve[ Complement[nlist, {n}] ][x /. n -> (-m), c] /;
    MemberQ[nlist, n] && FreeQ[m, n]
Convolve/: Convolve[nlist_List] [c_. Impulse[n_ + m_.], x_] :=
    Convolve[ Complement[nlist, {n}] ][c, x /. n -> (-m)] /;
    MemberQ[nlist, n] && FreeQ[m, n]

Convolve/: Convolve[n_][x_, z_?ZeroQ] := z
Convolve/: Convolve[n_][z_?ZeroQ, x_] := z

    (*  Formatting the convolution operator *)

Convolve/: Format[Convolve[n_]] :=
        Format[ Subscripted[Convolve[n]] ]
Convolve/: Format[Convolve[n_][x1_, x2_], TeXForm] :=
        StringForm[ "`` \\star_{``} ``", Format[x1], n, Format[x2] ]
Convolve/: Format[Convolve[n_][x1_, x2_, rest__], TeXForm] :=
        Block [    {format, i, length, operstr, xlist},
            xlist = {x1, x2, rest};
            length = Length[xlist];
            operstr = ToString[StringForm["\\star_{``}", n]];
            str = StringJoin[" ", operstr, " ``"];
            format = Apply[ StringJoin,
                    Table[ str, {i, 1, length} ] ];
            format = StringJoin["`` ", operstr, " ``", format];
            Apply[StringForm, {format, Map[Format, xlist]}] ]

(*  Delta  *)
Delta/: Delta[Infinity] := 0                (* definition         *)
Delta/: Delta[t_] := 0        /; N[ t != 0 ]
Delta/: Delta[-Infinity] := 0

Delta/: Delta[a_. x_Symbol + b_.] Delta[c_. x_ + d_.] := 0 /; ( b c != a d )
Delta/: Delta[a_. x_Symbol + b_.] Delta[a_. x_ + d_.] := 0 /; ( b != d )

Delta/: Extent1D[ Delta[a_. t_ + b_.], t_ ] := { -b/a, -b/a } /;
    FreeQ[{a,b}, t]

Delta/: Delta[t_, args__] := Delta[t] Delta[args]    (* multidimensional  *)
Delta/: Delta[List[args__]] := Delta[args]        (* separability         *)

Delta/: Format[ Delta, TeXForm ] := "\\delta"        (* output form       *)

Derivative[1][Delta] := Unit[1]                (* derivative rule   *)

    (* The following Delta rules were written by John Mitchell *)
    (* and modified by Brian Evans.                   *)
getivars[ x_Symbol ] := x
getivars[ {x_Symbol, xl_, xu_} ] := x
Unprotect[Integrate]                    (* integration rules *)
Integrate/: Integrate[a_. + expr_. Delta[term_], vars_] :=
    Integrate[a, vars] + DeltaIntegrate[expr Delta[term], vars] /;
    ! FreeQ[ term, getivars[vars] ]
Integrate/: Integrate[a_. + expr_. Delta[term_], vars1___, vars_, vars2___] :=
    Integrate[a, vars1, vars, vars2] +
      Integrate[ DeltaIntegrate[expr Delta[term], vars], vars1, vars2 ] /;
    ! FreeQ[ term, getivars[vars] ]
Protect[Integrate]

(*  DeltaIntegrate  *)
DeltaIntegrate[expr_. Delta[a_. x_ + b_.], x_Symbol] :=
    Limit[expr, x -> -b/a] CStep[x + b/a] / Abs[a] /;
    FreeQ[{a,b}, x] && FreeQ[expr, Delta]
DeltaIntegrate[expr_. Delta[a_. x_ + b_.], {x_Symbol, x1_, x1_}] := 0
DeltaIntegrate[expr_. Delta[a_ x_ + b_.], {x_Symbol, x1_, x2_}] :=
    If [ N[x1 > x2],
         -1 DeltaIntegrate[expr Delta[a x + b], {x, x2, x1}],
         Integrate[expr Delta[x + b/a] / Abs[a], {x, x1, x2}],
         Integrate[expr Delta[x + b/a] / Abs[a], {x, x1, x2}] ] /;
    FreeQ[{a,b}, x] && FreeQ[expr, Delta]
DeltaIntegrate[expr_. Delta[x_ + c_.], {x_Symbol, x1_, x2_}] :=
    If [ N[x1 > x2],
         -1 DeltaIntegrate[expr Delta[x + c], {x, x2, x1}],
         Limit[expr, x -> -c] CStep[-c - x1] CStep[c + x2],
         Limit[expr, x -> -c] CStep[-c - x1] CStep[c + x2] ] /;
    FreeQ[c, x] && FreeQ[expr, Delta]
        (* CStep[-c - x1] CStep[c + x2] <--> x1 <= -c <= x2 *)
DeltaIntegrate[Delta[a_. x_ + b_.] Delta[c_.x_ + d_.],
              {x_Symbol, x1_, x2_}] :=
    If [ N[x1 > x2],
         -1 DeltaIntegrate[Delta[a x + b] Delta[c x + d], {x, x2, x1}],
         (1/Abs[a c]) Delta[-b/a + d/c] *
        CStep[-(b/a) - x1] CStep[(b/a) + x2],
         (1/Abs[a c]) Delta[-b/a + d/c] *
        CStep[-(b/a) - x1] CStep[(b/a) + x2] ] /;
    FreeQ[{a,b,c,d}, x]
DeltaIntegrate[Delta[a_. x_ + b_.]^2, {x_Symbol, x1_, x2_}] :=
    If [ N[x1 > x2],
         -1 DeltaIntegrate[Delta[a x + b]^2, {x, x2, x1}],
         (1/a^2) Delta[0] CStep[-(b/a) - x1] CStep[(b/a) + x2],
         (1/a^2) Delta[0] CStep[-(b/a) - x1] CStep[(b/a) + x2] ] /;
    FreeQ[{a,b}, x]

(*  DeltaPlot--  plot Delta functions as special case of continuous plotting  *)

rewriteSymbolicSummation[index_, istart_, iend_, inc_,
             c_, b_, a_, {t_, tmin_, tmax_} ] :=
        Block [ {aa = a, adjfun, bb = b, cc = c, defaults, f1, f2,
           i, indexfunoft, minindex, maxindex, rules, varlist},
          varlist = Select[GetVariables[{a, b, c}],
                      ((! SameQ[#1, t]) &&
                    (! SameQ[#1, index]))&];
          If [ ! EmptyQ[varlist],
               defaults = Map[(#1 -> 1)&, varlist];
               aa = a /. defaults;
               bb = b /. defaults;
               cc = c /. defaults ];
          rules = Solve[ bb t + aa == 0, index ];
          indexfunoft = Replace[index, First[rules] ];
          f1 = N[ indexfunoft /. t -> tmin ];
          f2 = N[ indexfunoft /. t -> tmax ];
          minindex = Max[ istart, Ceiling[Min[f1, f2]] ];
          maxindex = Min[ iend, Max[f1, f2] ];
          adjfun = ((cc /. t -> - aa / bb) Delta[bb t + aa]) /.
               index -> i;
          Apply[ Plus,
             Table[adjfun, {i, minindex, maxindex, inc}] ] ]

dpcleanup[deltafuns_, t_, tmin_, tmax_] :=
    Block [    {delfuns, itemp, varlist = {}},
        delfuns = deltafuns /.
                ( d_ Summation[index_, istart_, iend_, inc_][c_] :>
                  Summation[index, istart, iend, inc][c d] );
        delfuns = delfuns /.
        { Summation[index_, istart_, iend_, inc_]
               [c_. Delta[b_. t + a_.]] :>
            rewriteSymbolicSummation[ index, istart, iend, inc, c,
                          b, a, {t, tmin, tmax} ],
          Periodic[k_, t][c_. Delta[b_. t + a_.]] :>
            rewriteSymbolicSummation[ itemp, -Infinity, Infinity, 1,
                          c, b, a + itemp k b,
                          {t, tmin, tmax} ] };
        varlist = Sort[ varlist ~Join~
                Select[ GetVariables[delfuns],
                    (! SameQ[#1, t])& ] ];
        If [ ! EmptyQ[varlist],
             Message[ DeltaPlot::unresolved, varlist ];
             delfuns = delfuns /. Map[(#1 -> 1)&, varlist] ];

        delfuns ]

getheight[ c_. Delta[a_. t_ + b_.], t_ ] := N[ c / Abs[a] ]

getlocation[ c_. Delta[a_. t_ + b_.], t_ ] := N[ - b / a ]

DeltaPlot[deltafuns_, t_, start_, end_] :=
    DeltaPlot[deltafuns, t, start, end, 0, 1]

DeltaPlot[deltafuns_, t_, start_, end_, plot_] := plot /;
    FreeQ[deltafuns, t]

    (* Use PlotRange as function to find range of plot [see MJ v1 #4 p47] *)
DeltaPlot[deltafuns_, t_, start_, end_, plot_] :=
    Block [    {xmax, xmin, ylist, ymax, ymin},
        {{xmin, xmax}, {ymin, ymax}} = PlotRange[plot];
        DeltaPlot[deltafuns, t, start, end, plot, ymin, ymax] ]

DeltaPlot[deltafuns_, t_, start_, end_, min_, max_] := NullPlot /;
    FreeQ[deltafuns, t]

    (* every version of DeltaPlot calls this one *)
DeltaPlot[deltafuns_, t_, start_, end_, min_, max_] :=
    Block [    {bottom, delfuns, height, heightList, length,
         locations, nend, nstart, width},
        nstart = N[start];
        nend = N[end];
        delfuns = Expand[ dpcleanup[deltafuns, t, nstart, nend] ];
        delfuns = If [ SameQ[Head[delfuns], Plus],
                   Apply[List, delfuns],
                   { delfuns } ];
        delfuns = Select[ delfuns,
                  N[LessEqual[nstart, getlocation[#, t], nend]]& ];
        height = N[max - min];
        width = nend - nstart;
        locations = Map[ getlocation[#1, t]&, delfuns];
        If [ min <= 0,
             bottom = 0; length = max,
             bottom = min; length = height ];
        length = Min[height, max];
        Which [ SameQ[locations, {}],
                  NullPlot,
            SameQ[$DeltaFunctionScaling, None],
                  Map[Arrow2D[{#1, bottom}, width, height, length]&,
                  locations ],
            True,
              heightList = Map[ getheight[#1, t]&, delfuns ];
              Map[ Arrow2D[{#[[1]],bottom}, width, height, #[[2]]]&,
                   Transpose[ {locations, heightList} ] ] ] ]

DeltaPlot[deltafuns_, t_, start_, end_, plot_, min_, max_] :=
    Block [    {deltaplot, ymax = max, ymin = min},
        If [ ymax == ymin,
             Which [ ymin > 0,  ymin = -ymax,
                 ymin == 0, ymax = 1,
                 ymin < 0,  ymax = -ymin ] ];
        deltaplot = DeltaPlot[deltafuns, t, start, end, ymin, ymax];
        If [ SameQ[deltaplot, NullPlot],
             plot,
             Prepend[ deltaplot, plot ] ] ]

(*  Difference  *)
Difference/: TheFunction[Difference[k_, n_][f_]] :=
    Block [    {jj, newf},
        newf = TheFunction[f] /. n ~ReplaceWith~ (n - jj);
        Apply[ Plus,
               Table[(-1)^jj Binomial[k,jj] newf, {jj, 0, k}] ] ] /;
    IntegerQ[k] && k > 0

Difference/: Extent1D[ Difference[k_, n_][f_], n_ ] := Extent1D[f, n]

Difference/: N[Difference[1, n_][f_]] := f - ( f /. n ~ReplaceWith~ (n - 1) )
Difference/: N[Difference[k_,n_][f_]] :=
    Apply[Plus,
          Table[(-1)^jj Binomial[k,jj] f /. n ~ReplaceWith~ (n - jj),
            {jj, 0, k}] ] /;
    IntegerQ[k] && k > 0

Difference/: Format[ Difference[k_, n_] ] := informat[Difference, k, n]

Difference/: Format[ Difference, TeXForm ] := "\\Delta"
Difference/: Format[ Difference[k_,n_], TeXForm ] :=
        inTeXformat[ Difference, k, n ]

Difference/: GetOperatorVariables[ Difference[k_, n_] ] := n

(*  Dirichlet  *)
Dirichlet/: Dirichlet[len_, 0] := len            (* definition    *)
Dirichlet/: Dirichlet[len_, k_?EvenQ Pi] := len (-1)^(len - 1)
Dirichlet/: Dirichlet[len_Integer, w_?NumberQ] :=
        Sin[len w / 2] / ( len Sin[w / 2] )

Dirichlet/: TheFunction[Dirichlet[len_, w_]] :=
        Sin[len w / 2] / ( len Sin[w / 2] )

Dirichlet/: N[ Dirichlet[len_, w_] ] := N [ TheFunction[Dirichlet[len, w]] ]

Dirichlet/: Format[ Dirichlet[l_, n_] ] :=        (* output forms     *)
        Format[ Dirichlet Subscript[l] [n] ]

(*  DiscreteGraphics  *)
oneline[xycoord_] :=
    ToCollection[ Line[{{xycoord[[1]], 0}, xycoord}], Point[xycoord] ]
DiscreteGraphics[coordlist_List, lthick_:0.008, pthick_:0.024 ] :=
    Graphics[ Join[ { Thickness[lthick], PointSize[pthick] },
            N [ Map[oneline, coordlist] ] ] ]

(*  Downsample [Myers, 219] [Bamberger, ch. 2]  *)

    (* Error checking of downsample parameters  *)

Downsample/: Downsample[m_List, n_List][f_] :=
    MyMessage[ Downsample::convert,
           Downsample[DiagonalMatrix[m], n][f],
           m ] /;
    VectorQ[m] && SameQ[Length[m], Length[n]] && ! PatternQ

Downsample/: Downsample[m_List, n_List][f_] :=
    MyMessage[Downsample::badmD, f] /;
    VectorQ[m] && ! SameQ[Length[m], Length[n]]

Downsample/: Downsample[m_List, n_List][f_] :=
    MyMessage[Downsample::badmD, f] /;
    MatrixQ[m] && ! Apply[SameQ, Append[Dimensions[m], Length[n]]]

Downsample/: Downsample[m_List, n_Symbol][f_] :=
    MyMessage[Downsample::badmD, f]

        (*  degenerate cases  *)

Downsample/: Downsample[1, n_Symbol][c_] := c
Downsample/: Downsample[-1, n_Symbol][c_] := Rev[n][c]

    (* convert upsample into a functional form  *)

Downsample/: TheFunction[Downsample[m_, n_Symbol][f_]] :=
    TheFunction[f] /. ( n -> m n )
Downsample/: TheFunction[Downsample[m_, n_List][f_]] :=
    TheFunction[f] /. ( n ~ReplaceWith~ (m . n) ) /;
    MatrixQ[m] && VectorQ[n]

Downsample/: N [ Downsample[m_, n_][f_] ] :=
    N [ TheFunction[Downsample[m,n][f]] ]

    (*  extent of downsample operation  *)

Downsample/: Extent1D[ Downsample[m_, n_][f_], n_ ] :=
    intdivide[ Extent1D[f, n], m ]

    (* formatting of the downsample operator *)

Downsample/: Format[ Downsample[m_, n_List] ] :=
    Block [    {bar},
        bar = TableForm[ Apply[ Table,
                    { "  |  ", Dimensions[n] } ] ];
        Format[ SequenceForm[ Downsample,
                      Subscript[ TableForm[m] ],
                      Subscript[ bar ],
                      Subscript[ TableForm[n] ] ] ] ]
Downsample/: Format[ Downsample[m_, n_] ] :=
    Format[ Subscripted[Downsample[m,n]] ]

Downsample/: Format[ Downsample, TeXForm ] := "\\downarrow"
Downsample/: Format[ Downsample[m_, n_], TeXForm ] :=
    inTeXformat[ Downsample, m, n ]

Downsample/: GetOperatorVariables[ Downsample[m_, n_] ] := n

(*  DFT  [Myers, 221]  *)
DFT/: DFT[L_, n_] := DFT[L, n, DummyVariables[Length[n], Global`k]]
DFT/: DFT[L_, n_, k_] [ InvDFT[L_, ki_, n_] [f_] ] :=
    If [ SameQ[k, ki], f, f /. ReplaceWith[ki, k] ] /;
    Length[k] == Length[ki]

DFT/: DFT^-1[a__] := InvDFT[a]

DFT/: Format[ DFT, TeXForm ] := "{\\cal DFT}"
DFT/: Format[ DFT[L_, n_, k_], TeXForm ] :=
    StringForm[ "{\\cal DFT}_{``, L = ``}", n, L]

DFT/: GetOperatorVariables[ DFT[L_, n_, k_] ] := n

DFT/: Extent1D[ DFT[L_, n_, k_][f_], n_ ] := {0, L}
DFT/: Extent1D[ DFT[L_, n_, k_][f_], k_ ] := {0, L}

(*  DTFT  *)
DTFT/: DTFT[n_] := DTFT[n, DummyVariables[Length[n], Global`w]]
DTFT/: DTFT[n_, w_] [ InvDTFT[wi_, n_] [f_] ] :=
    If [ SameQ[w, wi], f, f /. ReplaceWith[wi, w] ] /;
    Length[w] == Length[wi]

DTFT/: DTFT^-1[a__] := InvDTFT[a]

DTFT/: Format[ DTFT, TeXForm ] := "{\\cal DTFT}"
DTFT/: Format[ DTFT[n_, w_], TeXForm ] := StringForm[ "{\\cal DTFT}_{``}", n ]

(*  Extent1D  *)
lastExpression = Null

modifiedQ[x_, n_] :=
    Block [    {},
        lastExpression = x /.
          ( f_ :> Collect[f, n] /;
            PolynomialQ[f, n] || MixedPolynomialQ[Expand[f], n] );
        ! SameQ[ x, lastExpression ] ]

Extent1DUnion[ {x1_, y1_}, {x2_, y2_} ] :=
    { mymin[x1, x2], mymax[y1, y2] }
Extent1DIntersection[ {x1_, y1_}, {x2_, y2_} ] :=
    { mymax[x1, x2], mymin[y1, y2] }
Extent1DOrder[ x1_, y1_ ] :=
    { mymin[x1, y1], mymax[x1, y1] }

Extent1D[t_][x_] := Extent1D[x, t]

Extent1D[x_ + y_, n_] := Extent1DUnion[ Extent1D[x, n], Extent1D[y, n] ]
Extent1D[x_ y_, n_] := Extent1DIntersection[ Extent1D[x, n], Extent1D[y, n] ]
Extent1D[0, n_] := {0, 0}
Extent1D[x_, n_] := Extent1D[lastExpression, n] /; modifiedQ[x, n]
Extent1D[x_, n_] := {-Infinity, Infinity}

(*  FT  *)
FT/: FT[n_] := FT[n, DummyVariables[Length[n], Global`w]]
FT/: FT[n_, w_] [ InvFT[wi_, n_] [f_] ] :=
    If [ SameQ[w, wi], f, f /. ReplaceWith[wi, w] ] /;
    Length[w] == Length[wi]

FT/: FT^-1[a__] := InvFT[a]

FT/: Format[ FT, TeXForm ] := "{\\cal F}"
FT/: Format[ FT[n_, w_], TeXForm ] := StringForm[ "{\\cal F}_{``}", n ]

(*  GetDeltaFunctions  *)
choose[{a_, aterm_}] := If [ SameQ[aterm, 0], a, aterm ]

    (* rules specific for 1-D functions of t  *)
GetDeltaFunctions[ c_. Delta[a_. t_Symbol + b_.], t_Symbol ] :=
    c Delta[a t + b]

    (* rules specific for m-D functions of t1, t2, ...  *)
GetDeltaFunctions[ c_. Delta[a_. t_Symbol + b_.], tlist_List ] :=
    c Delta[a t + b] /;
    MemberQ[tlist, t]

    (* rules for all dimensions  *)
GetDeltaFunctions[ a_ + b_, rest__ ] :=
    GetDeltaFunctions[a, rest] + GetDeltaFunctions[b, rest]
GetDeltaFunctions[ h_[b__], rest__ ] :=
    Block [ {arglist, deltalist},
        arglist = {b};
        deltalist = Map[ GetDeltaFunctions[#1, rest]&, arglist ];
        If [ SameQ[Apply[Plus, deltalist], 0],
             0,
             Apply[h, Map[choose, Transpose[{arglist, deltalist}]]] ] ]
GetDeltaFunctions[ x_, rest__ ] := 0


(*  Impulse  *)
Impulse/: Impulse[Infinity] := 0               (* definition       *)
Impulse/: Impulse[n_] := 0    /; N[ n != 0 ]
Impulse/: Impulse[n_] := 1    /; N[ n == 0 ]
Impulse/: Impulse[-Infinity] := 0

Impulse/: Extent1D[ Impulse[a_. n_ + b_.], n_ ] := { -b/a, -b/a } /;
    FreeQ[{a,b}, n]

Impulse/: Impulse[n_, args__] := Impulse[n] Impulse[args] (* multidimensional *)
Impulse/: Impulse[List[args__]] := Impulse[args]      (* separability     *)

Impulse/: Format[ Impulse[n_], TeXForm ] :=          (* output form rule *)
      StringForm["\\delta[``]", n]

(*  Interleave [Myers, 219-220]  *)
Interleave/: Interleave[n_][x_] := x               (* definition          *)
Interleave/: TheFunction[Interleave[n_][x1_, x__]] :=
    Block [    {index = Unique["k"], L, xlist},
        xlist = List[x1, x];
        L = Length[xlist];
        Summation[index, 0, L-1, 1]
             [ Shift[index, n] [Upsample[L, n]
                         [ Element[xlist, index + 1] ] ] ] ]
Interleave/: N [ Interleave[n_][x1_, x__] ] :=
    N [ TheFunction[Interleave[n][x1,x]] ]

Interleave/: Format[ Interleave[n_] ] :=        (* output form rules *)
        Format[ Subscripted[Interleave[n]] ]
Interleave/: Format[ Interleave[n_], TeXForm ] :=
        Format[ Subscripted[Interleave[n]] ]


(*  IntervalsToFunction  *)
IntervalsToFunction[interlist_, n_:Global`n, step_:Step, pulse_:Pulse] :=
    Block [    {int, len, signal},

        MakeOne[interval_] :=
          Block [
            {f, nminus, nplus, pulselen},
            f = interval[[1]];
            nminus = interval[[2]];
            nplus = interval[[3]];
            pulselen = nplus - nminus + 1;
            Which [ SameQ[nminus, -Infinity],
                  f step[nplus - n],
                SameQ[nplus, Infinity],
                  f step[n - nminus],
                True,
                  f pulse[pulselen, n - nminus] ] ];

        len = Length[interlist];
        signal = MakeOne[ interlist[[1]] ];
        For [ int = 2, int <= len, int++,
              signal += MakeOne[ interlist[[int]] ] ];

        signal ]

(*  InvDFT  [Myers, 221]  *)
InvDFT/: InvDFT[L_, k_] :=
    InvDFT[L, k, DummyVariables[Length[k], Global`n]]
InvDFT/: InvDFT[L_, k_, ni_] [ DFT[L_, n_, k_] [f_] ] :=
    If [ SameQ[ni, n], f, f /. ReplaceWith[n, ni] ] /;
    Length[ni] == Length[n]

InvDFT/: Format[ InvDFT, TeXForm ] := "{\\cal DFT}^{-1}"
InvDFT/: Format[ InvDFT[L_, k_, n_], TeXForm ] :=
      StringForm[ "{\\cal DFT}^{-1}_{``, L = ``}", k, L ]

InvDFT/: GetOperatorVariables[ InvDFT[L_, k_, n_] ] := k

InvDFT/: Extent1D[ InvDFT[L_, k_, n_][f_], n_ ] := {0, L}
InvDFT/: Extent1D[ InvDFT[L_, k_, n_][f_], k_ ] := {0, L}

(*  InvDTFT  *)
InvDTFT/: InvDTFT[w_] := InvDTFT[w, DummyVariables[Length[w], Global`n]]
InvDTFT/: InvDTFT[w_, ni_] [ DTFT[n_, w_] [f_] ] :=
    If [ SameQ[ni, n], f, f /. ReplaceWith[n, ni] ] /;
    Length[ni] == Length[n]

InvDTFT/: Format[ InvDTFT, TeXForm ] := "{\\cal DTFT}^{-1}"
InvDTFT/: Format[ InvDTFT[w_, n_], TeXForm ] :=
      StringForm[ "{\\cal DTFT}^{-1}_{``}", w ]

(*  InvFT  *)
InvFT/: InvFT[w_] := InvFT[w, DummyVariables[Length[w], Global`t]]
InvFT/: InvFT[w_, ti_] [ FT[t_, w_] [f_] ] :=
    If [ SameQ[ti, t], f, f /. ReplaceWith[t, ti] ] /;
    Length[ti] == Length[t]

InvFT/: Format[ InvFT, TeXForm ] := "{\\cal F}^{-1}"
InvFT/: Format[ InvFT[w_, t_], TeXForm ] :=
      StringForm[ "{\\cal F}^{-1}_{``}", w ]

(*  InvL  *)
InvL/: InvL[s_] := InvL[s, DummyVariables[Length[s], Global`t]]
InvL/: InvL[s_, ti_] [ L[t_, s_] [f_] ] :=
    If [ SameQ[ti, t], f, f /. ReplaceWith[t, ti] ] /;
    Length[t] == Length[ti]

InvL/: Format[ InvL, TeXForm ] := "{\\cal L}^{-1}"
InvL/: Format[ InvL[s_, t_], TeXForm ] := StringForm[ "{\\cal L}^{-1}_{``}", s ]

(*  InvZ  [Myers, 221]  *)
InvZ/: InvZ[zi_] := InvZ[zi, DummyVariables[Length[zi], Global`n]]
InvZ/: InvZ[z_, ni_] [ Z[n_, z_] [f_] ] :=
    If [ SameQ[ni, n], f, f /. ReplaceWith[n, ni] ] /;
    Length[n] == Length[ni]

InvZ/: Format[ InvZ[z_, n_], TeXForm ] := "{\\cal Z}^{-1}"
InvZ/: Format[ InvZ[z_, n_], TeXForm ] := StringForm[ "{\\cal Z}^{-1}_{``}", z ]

(*  L  *)
L/: L[t_] := L[t, DummyVariables[Length[t], Global`s]]
L/: L[t_, s_] [ InvL[si_, t_] [f_] ] :=
    If [ SameQ[s, si], f, f /. ReplaceWith[si, s] ] /;
    Length[s] == Length[si]

L/: L^-1[a__] := InvL[a]

L/: Format[ L, TeXForm ] := "{\\cal L}"
L/: Format[ L[t_, s_], TeXForm ] := StringForm[ "{\\cal L}_{``}", t ]

(*  LineImpulse  *)
Unprotect[samecoefficient];
Unprotect[mergecoefficient];

SetAttributes[samecoefficient, Listable];
samecoefficient[c1_, c2_] :=
    If [ SameQ[c1,c2],
         ! SameQ[c1, Null],
         SameQ[c1, Null] || SameQ[c2, Null] ]

SetAttributes[mergecoefficient, Listable];
mergecoefficient[c1_, c2_] :=
    Which [ SameQ[c1, Null], c2,
        SameQ[c2, Null], c1,
        True,            c1 ]

Protect[samecoefficient];
Protect[mergecoefficient];

LineImpulseMerge[nlist_, ncoeffs_, klist_, kcoeffs_] :=
    Block [    {commonvars, kscale, nscale, scaledkcoeffs, scaledncoeffs},
        scaledncoeffs = ncoeffs;
        scaledkcoeffs = kcoeffs;
        commonvars = Intersection[nlist, klist];
        If [ ! EmptyQ[commonvars],
             normvar = First[commonvars];
             nscale = AssociateItem[normvar, nlist, ncoeffs];
             scaledncoeffs = ncoeffs / nscale;
             kscale = AssociateItem[normvar, klist, kcoeffs];
             scaledkcoeffs = kcoeffs / kscale ];
        allvars = Union[nlist, klist];
        coeffs1 = AssociateItem[allvars, nlist, scaledncoeffs];
        coeffs2 = AssociateItem[allvars, klist, scaledkcoeffs];
        { Apply[And, samecoefficient[coeffs1, coeffs2]],
          LineImpulse[allvars, mergecoefficient[coeffs1, coeffs2]] } ]

LineImpulse/: LineImpulse[n_?AtomQ, c_] := Impulse[c n]

LineImpulse/: LineImpulse[nlist_, ncoeffs_] LineImpulse[klist_, kcoeffs_] :=
    LineImpulseMerge[nlist, ncoeffs, klist, kcoeffs] [[2]] /;
    ! EmptyQ[Intersection[nlist,klist]] && LineImpulseMerge[nlist, ncoeffs, klist, kcoeffs] [[1]]

LineImpulse/: TheFunction[ LineImpulse[nlist_, ncoeffs_] ] :=
    Impulse[ Apply[Plus, nlist, ncoeffs] ]

LineImpulse/: N [ LineImpulse[nlist_, ncoeffs_] ] :=
    N [ TheFunction[ LineImpulse[nlist, ncoeffs] ] ]

LineImpulse/: Format[ LineImpulse[nlist_List, ncoeffs_List] ] :=
    Format[Impulse[Apply[Plus, nlist ncoeffs]]] /;
    Length[nlist] == Length[ncoeffs]

(*  OperatorInOperExpr  *)
OperatorInOperExpr[ h_[a__][s__] ] := h[a]
OperatorInOperExpr[ h_[s__] ] := h

(*  OperatorName  *)
OperatorName[ h_[a__][s__] ] := h
OperatorName[ h_[a__] ] := h
OperatorName[ x_ ] := x

(*  ParametersInOperExpr  *)
ParametersInOperExpr[ h_[a__][s__] ] := a
ParametersInOperExpr[ h_[s__] ] := Null

(*  Periodic  *)
Periodic/: Format[ Periodic[k_, n_] ] := informat[Periodic, k, n]
Periodic/: Format[ Periodic[k_, n_], TeXForm ] := inTeXformat[Perodic, k, n]

Periodic/: GetOperatorVariables[ Periodic[k_, n_] ] := n

(*  PolyphaseDownsample  *)
PolyphaseDownsample/: GetOperatorVariables[ PolyphaseDownsample[m_, h_, n_] ] := n

(*  PolyphaseUpsample  *)
PolyphaseUpsample/: GetOperatorVariables[ PolyphaseUpsample[l_, h_, n_] ] := n

(*  Pulse  *)
Pulse[l_] := Pulse[l, Global`n]                (* default args    *)

Pulse[0, n_] := 0                     (* definition      *)
Pulse[L_, 0] := 1
Pulse[1, n_] := Impulse[n]
Pulse[Infinity, n_] := Step[n]
Pulse/: Pulse[l_, n_] := 1 /; N[0 <= n < l]
Pulse/: Pulse[l_, n_] := 0 /; N[( n < 0 ) || ( n >= l )]

Pulse/: TheFunction[Pulse[l_, n_]] := Step[n] - Step[n - l]

Pulse/: Extent1D[ Pulse[l_, a_. n_ + b_.], n_ ] :=
        Extent1DOrder[-b/a, (l-b)/a] /;
        FreeQ[{a,b,l}, n]

Pulse/: Format[ Pulse[l_, n_] ] :=            (* output forms    *)
        Format[ Pulse Subscript[l] [n] ]

Pulse/: Format[ Pulse, TeXForm ] := "\\sqcap"
Pulse/: Format[ Pulse[l_, n_], TeXForm ] :=
        StringForm[ "\\sqcap_{``} [``]", l, n ]

Derivative[0,1][Pulse] :=                (* derivative rule *)
    ( Impulse[#2] + Impulse[#2 - #1] )&

(*  RationalGCD  *)
RationalGCD[ ] := 0
RationalGCD[ x_ ] := rationalGCD[{x}]
RationalGCD[ arg_, rest__ ] := rationalGCD[ {arg, rest} ]

rationalGCD[ numlist_List ] :=
    Block [    {denomgcd, gcdfun, gcdlist, newlist, numergcd},
        newlist = Flatten[numlist];
        gcdlist = Map[ gcdfilter, newlist ];
        gcdfun = If [ TrueQ[$VersionNumber >= 2.0],
                  PolynomialGCD, GCD ];
        numergcd = Apply[gcdfun, Map[Numerator, gcdlist]];
        denomgcd = Apply[gcdfun, Map[Denominator, gcdlist]];
        numergcd / denomgcd ]

(*  Rev  *)
Rev/: Rev[n_][Rev[n_][x_]] := x
Rev/: Rev[n_][Rev[m_][x_]] :=
    Block [    {remaining},
        remaining = SetExclusion[ToList[n], ToList[m]];
        Which [ EmptyQ[remaining],
              x,
            SameQ[Length[remaining], 1],
              Rev[ remaining[[1]] ] [x],
            True,
              Rev[remaining][x] ] ] /;
    ! EmptyQ[Intersection[ToList[n], ToList[m]]]

Rev/: Extent1D[ Rev[n_][f_], n_ ] :=
    Block [    {extent},
        extent = Extent1D[f, n];
        { -Second[extent], -First[extent] } ]

Rev/: TheFunction[Rev[n_][x_]] := ( TheFunction[x] /. n ~ReplaceWith~ (-n) )
Rev/: N[ Rev[n_][x_] ] := N [ TheFunction[Rev[n][x]] ]

Rev/: Format[ Rev[n_] ] := Format[ Subscripted[Rev[n]] ]
Rev/: Format[ Rev[n_], TeXForm ] := Format[ Subscripted[Rev[n]] ]

(*  RewriteSummations  *)
donotrewrite[f_, t_] := FreeQ[f, Summation] && FreeQ[f, Periodic[k_,t] ]

rewriteFindLimits[ x_, t_, tlower_, tupper_ ] :=
    Block [    {ext, tl, tu},
        ext = Extent1D[x, t];
        tl = mymax[ext[[1]], tlower];
        tu = mymin[ext[[2]], tupper];
        {tl, tu} ]

RewriteSummations[ f_, t_, tl_, tu_ ] := f /; donotrewrite[f, t]

RewriteSummations[ f_, t_, tlower_, tupper_ ] :=
    Block [    {ext, ii, tl, tu},
        SPSimplify[
          f /.
            { x_. Periodic[k_, v_][g_] :>         (* use v_ not t *)
            Block [ {},
              {tl, tu} = rewriteFindLimits[x, t, tlower, tupper];
              x rewriteonesummation[ Simplify[g /. (t -> t + ii k)],
                             True, t, tl, tu, ii,
                         -Infinity, Infinity, 1 ] ] /;
             SameQ[v, t] && donotrewrite[g, t],
              x_. Summation[i_, il_, iu_, inc_][g_] :>
            Block [ {},
              {tl, tu} = rewriteFindLimits[x, t, tlower, tupper];
              x rewriteonesummation[ Simplify[g], False, t, tl, tu,
                         i, il, iu, inc ] ] /;
             donotrewrite[g, t] },
          Variables -> t ] ]

rewriteonesummation[ a_ + b_, flag_, t_, tl_, tu_, i_, il_, iu_, inc_ ] :=
    rewriteonesummation[a, flag, t, tl, tu, i, il, iu, inc] +
    rewriteonesummation[b, flag, t, tl, tu, i, il, iu, inc]

rewriteonesummation[ g_, periodicflag_, t_, tl_, tu_, i_, il_, iu_, inc_ ] :=
    Block [ {ext, il0, ilact, intflag, iu0, iuact,
         lext, lsol, retval, solutions, uext, usol, valid},
        intflag = IntegerQ[inc] &&    (* index var i is an integer *)
              ( IntegerQ[il] || SameQ[il, -Infinity] ) &&
              ( IntegerQ[iu] || SameQ[iu, Infinity] );
        ext = Simplify[ Extent1D[g, t] ] // minmaxSimplify;
        valid = ListQ[ext] && SameQ[Length[ext], 2];
        If [ ! valid, ext = { -Infinity, Infinity } ];

        {lext, uext} = ext;
        If [ FreeQ[lext, i] && ! FreeQ[uext, i], lext = uext ];
        If [ FreeQ[uext, i] && ! FreeQ[lext, i], uext = lext ];

        If [ FreeQ[lext, i] && FreeQ[uext, i],
             Message[ Summation::extent ];
               retval = If [ periodicflag,
                     g /. i -> 0,
                     Summation[i, il, iu, inc][g] ],

             lsol = Solve[ lext == tu, i ];
             usol = Solve[ uext == tl, i ];
             Which [ SameQ[Head[lsol], Solve] ||
                 SameQ[Head[usol], Solve],
                   Message[ Summation::extent ];
                     retval = If [ periodicflag,
                          g /. i -> 0,
                          Summation[i, il, iu, inc][g] ],

                 intflag,
                   solutions = N[ Join[ i /. lsol, i /. usol ] ];
                       il0 = Min[solutions];
                       iu0 = Max[solutions];
                       ilact = Max[il, Ceiling[Min[il0, iu0]]];
                       iuact = Min[iu, Max[il0, iu0]];
                       retval = Sum[ g, {i, ilact, iuact, inc} ],

                 False,
                   solutions = Join[ i /. lsol, i /. usol ];
                       il0 = Apply[Min, solutions] // minmaxSimplify;
                       iu0 = Apply[Max, solutions] // minmaxSimplify;
                       ilact = Max[il, Min[il0, iu0]] // minmaxSimplify;
                       iuact = Min[iu, Max[il0, iu0]] // minmaxSimplify;
                   If [ InfinityQ[il],
                            ilact = ilact - Mod[ilact, inc];
                              iuact = iuact - Mod[iuact, inc],
                            ilact = ilact - Mod[ilact - il, inc];
                              iuact = iuact - Mod[iuact - il, inc] ];
                       retval = Sum[ g, {i, ilact, iuact, inc} ] ] ];
    retval ]


(*  ScaleAxis [Myers, 220]  *)
ScaleAxis/: TheFunction[ScaleAxis[l_, w_Symbol][x_]] :=
        TheFunction[x] /. w -> l w
ScaleAxis/: TheFunction[ScaleAxis[{l__},{w__}][x_]] :=
        ( TheFunction[x] /. {w} ~ReplaceWith~ ({l} . {w}) )

ScaleAxis/: N[ ScaleAxis[l_, w_][x_] ] := N[ TheFunction[ ScaleAxis[l,w][x] ] ]

ScaleAxis/: Format[ ScaleAxis[l_,w_] ] := informat[ScaleAxis, l, w]
ScaleAxis/: Format[ ScaleAxis[l_,w_], TeXForm ] := inTeXformat[ScaleAxis, l, w]

ScaleAxis/: GetOperatorVariables[ ScaleAxis[l_, w_] ] := w

(*  ScalingFactor --  must be independent of variable z  *)
ScalingFactor[f_, z_] :=
    breakup[ Apply[ RationalGCD, GetAllFactors[f, z] ], z ]

(*  SequencePlot  *)
Options[ SequencePlot ] :=
    { Axes -> axesdefaultvalue,
      DisplayFunction :> $DisplayFunction,
      PlotLabel -> "Sequence Plot",
      PlotRange -> All }

baseticks = { 1, 2, 2, 5, 5, 5, 5, 5, 10, 10 };
getticks[nstart_, nend_] :=
    Block [    {factor, i, ninc, ns, nterm},
        ninc = Round[(nend - nstart) / 5];
        If [ ninc == 0, ninc = 1 ];
        If [ ninc > 2,
             factor = 10;
             While [ ninc > factor, factor *= 10 ];
             factor /= 10;
             nterm = intdivide[ninc, factor];
             ninc = factor * baseticks[[nterm]] ];
        ns = Sign[nstart] ninc intdivide[Abs[nstart], ninc];
        Table[i, {i, ns, nend, ninc}] ]

dgetvalue[f_, n_Symbol, n0_Integer] := GetValue[f, n, n0]
dgetvalue[f_, n_Symbol, n0_Real] := GetValue[f, n, Round[n0]]
dgetvalue[f_, {n1_Symbol, n2_Symbol}, {n01_?NumberQ, n02_?NumberQ}] :=
    GetValue[f, {n1, n2}, {Round[n01], Round[n02]}]

SequencePlot[f_List, {n1_Symbol, st1_, end1_},
             {n2_Symbol, st2_, end2_}, op___] :=
    MyMessage[ SequencePlot::nolist, NullPlot ]

SequencePlot[f_, {n1_Symbol, st1_, end1_}, {n2_Symbol, st2_, end2_}, op___] :=
    Block [    {funct, ftemp, funs, n1length, n2length, n1str, n2str,
         n1temp, n2temp, nstart, nend, oplist, plot},

        nstart = Floor[N[{st1, st2}]];
        nend   = Ceiling[N[{end1, end2}]];
        n1length = nend[[1]] - nstart[[1]];
        n2length = nend[[2]] - nstart[[2]];
        n1str = StripPackage[n1];
        n2str = StripPackage[n2];
        oplist = Join [    ToList[op],
                Options[SequencePlot],
                { PlotPoints -> {n1length + 1, n2length + 1},
                  AxesLabel -> {n1str, n2str, " "},
                  Ticks -> Automatic } ];

        funct = ToDiscrete[f];
        ftemp = N[ TheFunction[ funct ] ];

        (*  Generate three-dimensional coordinates for plot  *)
        Off[Power::infy, Infinity::indt];
        plot = Apply[ Plot3D,
                  Join[ { dgetvalue[ ftemp, {n1, n2},
                         {n1temp, n2temp} ],
                      { n1temp, nstart[[1]], nend[[1]] },
                      { n2temp, nstart[[2]], nend[[2]] } },
                    oplist ] ];
        On[Power::infy, Infinity::indt];

        plot ]


SequencePlot[f_, {n_Symbol, start_, end_}, op___] :=
    Block [    {factor, ftemp, funct, funs, length, linethickness,
         nend, nstart, nstr, ntemp, oplist, plot,
         pointthickness, showoptions},

        nstart = Floor[N[start]];
        nend   = Ceiling[N[end]];

        (*  Set up the options  *)
        nticks = getticks[nstart, nend];
        nstr = StripPackage[n];
        oplist = Join [    Flatten[ToList[op]],
                Options[SequencePlot],
                    { AxesLabel -> { nstr, " " },
                  Ticks -> { nticks, Automatic } } ];

        factor = N[ Log[10, nend - nstart + 3 ] ];
        linethickness = 0.008 / factor;
        pointthickness = 0.024 / factor;

        (*  Extract the functional form of the expression  *)
        funct = ToDiscrete[f];
        ftemp = N[ TheFunction[ funct ] ];

        (*  Handle any infinite summations  *)
        ftemp = RewriteSummations[ftemp, n, nstart, nend];

        (*  Generate form for 1-D function  *)
        SetAttributes[seqFun1D, Listable];    (* avoid fun pointers *)
        seqFun1D[g_] :=
            DiscreteGraphics [
                Table [ {ntemp, dgetvalue[g, n, ntemp]},
                    {ntemp, nstart, nend}],
                linethickness,
                pointthickness ];

        (*  Plot "time"-domain    *)
        Off[Power::infy, Infinity::indt];
        showoptions = RemoveOptions[oplist, specialPlotOptions];
        plot = Show [ seqFun1D[ftemp], showoptions ];
        On[Power::infy, Infinity::indt];

        plot ]

(*  SequenceToFunction  *)
SequenceToFunction[seq_, n_:Global`n, noffset_:0] :=
    Block [ {},
        (* Variable no. of arguments -- could not use Function *)
        impulsegen[l__] := seq[[l]] Apply[Impulse,
                          n - {l} - noffset + 1];
        Apply[ Plus, Flatten [ Array[impulsegen, Dimensions[seq]] ] ] ]

(*  SignalCleanup  *)
SignalCleanup[funct_, t_Symbol, nstart_, nend_] :=
    Block [    {deltafuns, ftemp},

        (*  Separate delta functions from rest of funct     *)
        (*  DeltaPlot handles Summation operators by itself *)

        deltafuns = GetDeltaFunctions[funct, t];
        ftemp = funct //.
            { h_[ c_. Delta[a_. t + b_.] ] :> 0,
              h_[ c_. Delta[a_. t + b_.] + x_ ] :> h[x],
              Delta[a_. t + b_.] :> 0 };

        (*  Convert function into plottable functional form  *)

        ftemp = N[ TheFunction[ftemp] ];

        (*  Rewrite infinite summations and periodic operators  *)

        deltafuns = RewriteSummations[deltafuns, t, nstart, nend];
        ftemp = RewriteSummations[ftemp, t, nstart, nend];

        { deltafuns, ftemp } ]


SignalCleanup[funct_, {t1_Symbol, t2_Symbol}, nstart_, nend_] :=
    Block [    {deltafuns, ftemp},

        (*  Separate delta functions from rest of funct  *)

        deltafuns = GetDeltaFunctions[funct, {t1, t2}];
        ftemp = funct //.
            { h_[ c_. Delta[a_. t1 + b_.] ] :> 0,
              h_[ c_. Delta[a_. t1 + b_.] + x_ ] :> h[x],
              Delta[a_. t1 + b_.] :> 0,
              h_[ c_. Delta[a_. t2 + b_.] ] :> 0,
              h_[ c_. Delta[a_. t2 + b_.] + x_ ] :> h[x],
              Delta[a_. t2 + b_.] :> 0 };

        (*  Convert freq. response into plottable function form  *)

        ftemp = N [ TheFunction[funct] ];

        { deltafuns, ftemp } ]

(*  SignalPlot  *)
specialPlotOptions =
    Complement[ Sort[ Map[First, Options[Plot]] ],
            Sort[ Map[First, Options[Graphics]] ] ]

sigPlot1D[args__][g_] := SignalPlot[g, args]

Options[ SignalPlot ] :=
    { Axes -> axesdefaultvalue,
      DisplayFunction :> $DisplayFunction,
      PlotLabel -> "Signal Plot",
      PlotRange -> All,
      Ticks -> Automatic }

SignalPlot[ f_List, {t1_Symbol, start1_, end1_},
            {t2_Symbol, start2_, end2_}, op___ ] :=
    MyMessage[ SignalPlot::nolist, NullPlot ]

SignalPlot[ f_,    {t1_Symbol, start1_, end1_},
        {t2_Symbol, start2_, end2_}, op___ ] :=
    Block [    {deltafuns, freqresp, ftemp, funct, oplist,
         t1str, t2str, t1temp, t2temp},

        (*  Set up for plotting  *)
        t1str = StripPackage[t1];
        t2str = StripPackage[t2];
        oplist = Join[ ToList[op],
                   { AxesLabel -> { t1str, t2str, " " } },
                   Options[ SignalPlot ] ];

        (*  Make function continuous and replace non-sep delta's *)
        funct = ToContinuous[f] /.
             Delta[a_. t1 + b_. t2 + c_.] :> Impulse[a t1 + b t2 + c];

        (*  Separate delta functions from rest of funct and      *)
        (* convert the non-delta expressions into plottable form *)
        { deltafuns, ftemp } =
            SignalCleanup[ funct, {t1,t2},
                       N[{start1, start2}], N[{end1, end2}] ];

        (*  Plot delta functions separately from rest of function  *)

        Plot3D[ GetValue[ftemp, {t1, t2}, {t1temp, t2temp}],
            { t1temp, start1, end1 },
            { t2temp, start2, end2 },
            Release[oplist] ] ]

SignalPlot[f_List, {t_Symbol, start_, end_}, op___] :=
    Block [    {fun, oplist, showoptions, theplots},
        fun = sigPlot1D[{t,start,end}, DisplayFunction->Identity, op];
        theplots = Map[ fun, f ];
        oplist = Join[ToList[op], Options[SignalPlot]];
        showoptions = RemoveOptions[oplist, specialPlotOptions];
        Show[ theplots, ToCollection[showoptions] ] ]

SignalPlot[f_, {t_Symbol, start_, end_}, op___] :=
    Block [    {deltafuns, ftemp, i, imstyle,
         nend = N[end], nstart = N[start], oplist,
         plot, showoptions, restyle, tstr, ttemp },

        (*  Set up for plotting  *)

        tstr = StripPackage[t];
        restyle = Thickness[0.006];
        imstyle = Dashing[{0.05,0.05}];
        oplist = Join[ ToList[op],
                   { AxesLabel -> { tstr, " " },
                 DisplayFunction -> Identity,
                 PlotStyle -> { restyle, imstyle } },
                   Options[ SignalPlot ] ];

        (*  Separate delta functions from rest of funct and  *)
        (*  convert function into plottable functional form  *)

        { deltafuns, ftemp } =
            SignalCleanup[ ToContinuous[f], t, nstart, end ];

        (*  Plot delta functions separately from ftemp  *)

        Off[Power::infy, Infinity::indt];

        (*  1.  Generate the graphics for the entire signal  *)

        If [ ZeroQ[ftemp],
             plot = If [ ! AtomQ[deltafuns],   (* delta funs found *)
                 DeltaPlot[deltafuns, t, nstart, nend],
                 NullPlot ];
             plot = { Plot[ 0,
                    { ttemp, nstart, nend},
                    Release[oplist] ],
                  plot },

             plot = Plot[ { Re[GetValue[ftemp, t, ttemp]],
                    Im[GetValue[ftemp, t, ttemp]] },
                  { ttemp, nstart, nend },
                  Release[oplist] ];
             If [ ! AtomQ[deltafuns],      (* delta funs found *)
                  plot = DeltaPlot[deltafuns, t, nstart, nend, plot]] ];

        (*  2.  Plot the entire signal  *)

        oplist = Join[ ToList[op], Options[SignalPlot] ];
        showoptions = RemoveOptions[oplist, specialPlotOptions];

        plot = Show[ plot, ToCollection[showoptions] ];

        On[Power::infy, Infinity::indt];

        plot ]


(*  SignalsInOperExpr  *)
SignalsInOperExpr[ h_[s__] ] := s
SignalsInOperExpr[ h_[a__][s__] ] := s

SignalsInOperExpr[ h_[s__], h_ ] := s
SignalsInOperExpr[ h_[s__], op_ ] := h[s]        /; ! SameQ[h, op]

SignalsInOperExpr[ h_[a__][s__], h_ ] := s
SignalsInOperExpr[ h_[a__][s__], op_ ] := h[a][s]    /; ! SameQ[h, op]

(*  Shift [Myers, 217]  *)
Shift/: Shift[m_, n_List] := Shift[ Table[m, Length[n]], n ] /; ! ListQ[m]

Shift/: Shift[m_List, n_List][x_] := MyMessage[ Shift::badmD, x ] /;
    Length[m] != Length[n]

Shift/: Shift[0, n_][x_] := x

Shift/: Extent1D[ Shift[m_, n_][f_], n_ ] := Extent1D[f, n] + m

Shift/:    TheFunction[Shift[k_,n_][x_]] :=
    ( TheFunction[x] /. n ~ReplaceWith~ (n - k) )
Shift/: N[ Shift[k_,n_][x_] ] := N [ TheFunction[ Shift[k,n][x] ] ]

Shift/: Format[ Shift[m_,n_] ] := informat[Shift, m, n]
Shift/: Format[ Shift[m_,n_], TeXForm ] := inTeXformat[Shift, m, n]

Shift/: GetOperatorVariables[ Shift[l_, n_] ] := n


(*  Sinc  *)
Sinc/: Sinc[-Infinity] := 0
Sinc/: Sinc[0] := 1
Sinc/: Sinc[0.] := 1.
Sinc/: Sinc[Infinity] := 0

Sinc/: Sinc[x_?NumberQ] := Sin[x] / x

Sinc/: Sinc[x_, y__] := Sinc[x] Sinc[y]

Sinc/: Format[Sinc, TeXForm] := "{\\rm sinc}"

Sinc/: Derivative[1][Sinc] := (( Cos[#1] / #1 ) - ( Sinc[#1] / #1 ))&


(*  SPSimplify  *)
Options[SPSimplify] :=
    Join[{ Dialogue -> False, Variables -> {} }, Options[Simplify]]

kludgeRules =   (* prevents infinite loop in Mathematica 2.0 and lower *)
      {    a1_. x_. Exp[Complex[0, i1_] th1_.] +
        a2_. x_. Exp[Complex[0, i2_] th2_.] :>
        2 x Abs[a1] Cos[Arg[a1] + i1 th1] /;
        ( SameQ[a1, Conjugate[a2]] && SameQ[i1 th1, -i2 th2] ),
        a1_. x_. Exp[Complex[0, i1_] th1_.] +
        a2_. x_. Exp[Complex[0, i2_] th2_.] :>
        2 I x Abs[a2] Sin[Arg[a2] + i2 th2] /;
        ( SameQ[a1, -a2] && SameQ[i1 th1, -i2 th2] ) }

SPSimplify[ x_, op___ ] :=
    Block [    {dialflag, different, laste, newe, oplist,
         trig, varrules, varflag, variables},
        oplist = ToList[op] ~Join~ Options[SPSimplify];
        trig = Replace[Trig, oplist];
        dialflag = InformUserQ[oplist];
        variables = Replace[Variables, oplist];
        varflag = ! EmptyQ[variables];
        If [ varflag,
             varrules = makesimplifyrules[variables] ];

        newe = x;
        different = True;
        While [ different,
            laste = newe;
            If [ ! EmptyQ[ kludgeRules ],
                 newe = newe /. kludgeRules ];
            newe = If [ TrueQ[ $VersionNumber >= 2.0 ],
                    Simplify[newe, Trig -> trig],
                    Simplify[newe] ];
            newe = newe /. SPSimplificationRules;
            If [ varflag, newe = newe /. varrules ];
            different = ! SameQ[laste, newe];
            If [ dialflag && different,
                 Print["which simplifies to"]; Print[newe] ] ];
        newe ]

makesimplifyrules[var_Symbol] := makesimplifyrules[var, {}]
makesimplifyrules[var1_Symbol, var2_Symbol, rest___] :=
    makesimplifyrules[{var1, var2, rest}]
makesimplifyrules[varlist_List] :=
    Apply[ Join,
           Map[ makesimplifyrules[#,
                      Complement[varlist, ToList[#]]]&,
                      varlist] ]
makesimplifyrules[t_, varlist_] :=
    { Delta[a_ t + b_.] :> Delta[t + b/a] / Abs[a] /;
        FreeQ[a, t] && MyFreeQ[b, varlist],
      d_. (y_. + amp1_. Delta[t + t0_.]) :>
        d y + Limit[TheFunction[d amp1], t -> -t0] Delta[t + t0] /;
        FreeQ[t0, t] && (! FreeQ[d amp1, t]) &&
        MyFreeQ[t0, varlist] &&
          ! InfinityQ[ Limit[TheFunction[d amp1], t -> -t0] ],
      c_. Sign[b_. t] + c_ :> 2 c CStep[b t] /;
        FreeQ[b, t],
      CStep[a_ t + b_.] :> CStep[t + b/a] / Abs[a] /;
        Positive[a] && FreeQ[a, t] && MyFreeQ[b, varlist],
      CStep[a_ t + b_.] :> CStep[-t - b/a] / Abs[a] /;
        Negative[a] && FreeQ[a, t] && MyFreeQ[b, varlist] }

(*  Step  *)
Step/: Step[Infinity] := 1                (* definition        *)
Step/: Step[t_] := 1        /; N[ t > 0 ]
Step/: Step[t_] := 0        /; N[ t < 0 ]
Step/: Step[0] := 1
Step/: Step[0.] := 1
Step/: Step[-Infinity] := 0

Step/: Extent1D[ Step[a_. n_ + b_.], n_ ] :=
        Extent1DOrder[-b/a, Sign[a] Infinity] /;
        FreeQ[{a,b}, n]

Step/: Step[t_, args__] := Step[t] Step[args]        (* multidimensional *)
Step/: Step[List[args__]] := Step[args]            (* separability        *)

Step/: Step[t_]^a_ := Step[t]                (* automatic        *)
Step/: Step[n_Symbol + m_.] Step[n_Symbol  + k_.] :=    (* simplification   *)
    Step[n + Min[m,k]] /;                (* rules        *)
    NumberQ[m] && NumberQ[k]
Step/: Step[n_Symbol + m_.] Step[-n_Symbol + k_.] :=
    Pulse[k + m, n + m] /;
    N[ -m <= k ]
Step/: Step[n_Symbol + m_.] Step[-n_Symbol + k_.] := 0 /;
    N[ -m > k ]

Step/: Simplify[Step[a_ n_Symbol + b_.]] :=        (* man. simp. rules *)
    Step[n + b/a] /;
    MyFreeQ[{a,b}, n] && Implies[ NumberQ[N[a]], N[a > 0] ]
Step/: Simplify[Step[- a_ n_Symbol + b_.]] :=
    Step[-n + b/a] /;
    MyFreeQ[{a,b}, n] && Implies[ NumberQ[N[a]], N[a > 0] ]

Step/: Format[ Step[x_], TeXForm ] :=            (* output form      *)
    StringForm["u[``]", x]

Derivative[1][Step] := Impulse                (* derivative rule  *)

(*  Summation  *)
Summation/: Summation[i_Symbol, n_] := Summation[i, 1, n, 1]
Summation/: Summation[i_Symbol, istart_, iend_] :=
    Summation[i, istart, iend, 1]
Summation/: Summation[i_Symbol, istart_, iend_, inc_][0] := 0

Summation/: TheFunction[Summation[i_Symbol, istart_, iend_, inc_][f_]] :=
    Block [    {tempsym},
        If [ NumberQ[istart] && NumberQ[iend] && NumberQ[inc],
             Sum[ TheFunction[f /. i -> tempsym],
              {tempsym, istart, iend, inc} ],
             Summation[i, istart, iend, inc][ TheFunction[f] ] ] ]

Summation/: Extent1D[ Summation[i_, istart_, iend_, inc_][f_], n_ ] :=
    Block [    {extent},
        extent = Extent1D[f, n];
        Extent1DOrder[ extent /. i -> istart,
                   extent /. i -> iend ] ] /;
    FreeQ[{i,istart,iend,inc}, n]

Summation/: N[Summation[i_Symbol, istart_?NumberQ, iend_?NumberQ, inc_?NumberQ][f_]] :=
    Sum[ N[f], {i, istart, iend, inc} ]

Summation/: Format[ Summation[i_, istart_, iend_, 1], TeXForm ] :=
    StringForm[ "\\sum_{`` = ``}^{``}", i, istart, iend ]

Summation/: Format[ Summation[i_, istart_, iend_, inc_], TeXForm ] :=
    StringForm["\\sum_{``=``, `` = `` + ``}^{``}",
           i, istart, i, i, inc, iend]

    (* The following rules were suggested by John Mitchell *)
getdvars[ x_Symbol ] := x
getdvars[ {x_Symbol, k_} ] := x
Unprotect[D]
D/: D[expr_. Summation[i_, istart_, iend_, inc_][f_], x__] :=
    Summation[i, istart, iend, inc][ D[expr f, x] ] /;
    MyFreeQ[ i, Map[getdvars, {x}] ]
Protect[D]
Unprotect[Integrate]
Integrate/: Integrate[ expr_. Summation[i_, istart_, iend_, inc_][f_], x__] :=
     Summation[i, istart, iend, inc][ Integrate[expr f, x] ] /;
     MyFreeQ[ i, Map[getdvars, {x}] ]
Protect[Integrate]

(*  TheFunction  *)
TheFunction[x_?AtomQ] := x
TheFunction[Hold[x__]] := Hold[x]
TheFunction[Sign[x_]] := CStep[x] - CStep[-x]    (* Step or CStep will work  *)
TheFunction[h_[a_]] := h [ TheFunction[a] ]
TheFunction[h_[a1_, a__]] := Apply[ h, Map[TheFunction, { a1, a }] ]

(*  ToContinuous  *)
ToContinuous[expr_] := ReplaceRepeated[expr, ToContinuousRules]

ToContinuousRules =
     {
    Impulse[n_]   :> Delta[n],
    Step[n_]      :> CStep[n],
    Pulse[l_, n_] :> CPulse[l, n],
    Convolve[n_]  :> CConvolve[n]
     }

(*  ToDiscrete  *)
ToDiscrete[expr_] := ReplaceRepeated[expr, ToDiscreteRules]

ToDiscreteRules =
     {
    CStep[t_]      :> Step[t],
    Delta[t_]      :> Impulse[t],
    CPulse[l_, t_] :> Pulse[l, t],
    CConvolve[t_]  :> Convolve[t]
     }

(*  Unit *)
Unit/: Unit[0][t_] := Delta[t]                (* definition       *)
Unit/: Unit[-1][t_] := CStep[t]
Unit/: Unit[n_][t_?NumberQ] :=
     t^-(n + 1) / Factorial[-(n + 1)] CStep[t] /; N[n < 0]

Unit/: Extent1D[ Unit[n_Integer][a_. t_ + b_.], t_ ] := {-b/a, -b/a} /;
    FreeQ[{a,b}, t] && ( n <= 0 )
Unit/: Extent1D[ Unit[n_Integer][a_. t_ + b_.], t_ ] :=
    Extent1DOrder[ -b/a, Sign[a] Infinity ] /;
    FreeQ[{a,b}, t] && ( n > 0 )

Unit/: Unit[n_][List[args]] := Unit[n][args]        (* multidimensional *)
Unit/: Unit[n_][t1_, t__] := Unit[n][t1] Unit[n][t]    (* separability     *)

Unit/: TheFunction[Unit[n_][t__]] :=
     t^-(n + 1) / Factorial[-(n + 1)] CStep[t] /; N[n < 0]

Unit/: Format[ Unit[n_] ] := Subscripted[ Unit[n] ]    (* output forms     *)
Unit/: Format[ Unit[n_], TeXForm ] := StringForm[ "u_{``}", n ]

Derivative[1][Unit[n_]] := Unit[n+1]            (* derivative rule  *)

Unprotect[Integrate]                    (* integration rules *)
Integrate/: Integrate[ expr_. Unit[n_][a_. x_ + b_.], x_Symbol ] :=
    Limit[expr, x -> -b/a] Unit[n-1][a x + b] / a /;
    FreeQ[{a,b,c}, x]
Integrate/: Integrate[ expr_. Unit[n_][a_ x_ + b_.], {x_Symbol, x1_, x2_}] :=
    Integrate[expr Unit[n][x + b/a] / a, x] /;
    FreeQ[{a,b}, x]
Integrate/: Integrate[expr_. Unit[n_][x_ + c_.], {x_Symbol, x1_, x2_}] :=
    Limit[expr, x -> -c] CStep[-c - x1] CStep[c + x2] Unit[n-1][x + c] /;
    FreeQ[c, x]    (* CStep[-c - x1] CStep[c + x2] <--> x1 <= -c <= x2 *)
Protect[Integrate]

(*  Upsample [Myers, 218-219] [Bamberger, ch. 2]  *)

    (* Error checking of upsample parameters  *)

Upsample/: Upsample[m_List, n_List][f_] :=
    MyMessage[Upsample::convert, Upsample[DiagonalMatrix[m], n][f], m] /;
    VectorQ[m] && SameQ[Length[m], Length[n]]

Upsample/: Upsample[m_List, n_List][f_] :=
    MyMessage[Upsample::badmD, f] /;
    VectorQ[m] && ! SameQ[Length[m], Length[n]]

Upsample/: Upsample[m_List, n_List][f_] :=
    MyMessage[Upsample::badmD, f] /;
    MatrixQ[m] && ! Apply[SameQ, Append[Dimensions[m], Length[n]]]

Upsample/: Upsample[m_List, n_Symbol][f_] :=
    MyMessage[Upsample::badmD, f]

        (*  degenerate cases  *)

Upsample/: Upsample[1, n_][f_] := f
Upsample/: Upsample[-1, n_][f_] := Rev[n][f]

    (* convert upsample into a functional form  *)

Upsample/: N[Upsample[m_,n_][f_List]] := UpsampleSequence[f, m, 0]

Upsample/: TheFunction[Upsample[m_, n_Symbol][f_List]] :=
        UpsampleSequence[f, m, 0]
Upsample/: TheFunction[Upsample[m_, n_Symbol][f_]] :=
        UpsampledFunction[m, 0, n, TheFunction[f] /. n -> n / m]
Upsample/: TheFunction[Upsample[m_, n_List][f_]] :=
        UpsampledFunction[ m, 0, n,
            TheFunction[f] /. ReplaceWith[n, Inverse[m] . n] ]

    (* extent of upsampled function  *)

Upsample/: Extent1D[ Upsample[m_, n_][f_], n_ ] := m Extent1D[f, n]

    (* formatting of the upsample operator *)

Upsample/: Format[ Upsample[m_, n_List] ] :=
    Block [    {bar},
        bar = TableForm[ Table[ "  |  ", Release[Dimensions[n]] ] ];
        Format[ SequenceForm[ Upsample,
                      Subscript[ TableForm[m] ],
                      Subscript[ bar ],
                      Subscript[ TableForm[n] ] ] ] ]
Upsample/: Format[ Upsample[m_, n_] ] :=
        Format[ Subscripted[Upsample[m, n]] ]

Upsample/: Format[ Upsample, TeXForm ] := "\\uparrow"
Upsample/: Format[ Upsample[m_, n_], TeXForm ] :=
        inTeXformat[ Upsample, m, n ]

Upsample/: GetOperatorVariables[ Upsample[l_, n_] ] := n

(*  UpsampledFunction  *)
lastMatrix = Null
lastMatrixInverse = Null

integerVectorQ[x_] := Apply[ And, Map[ IntegerQ, x ] ]

UpsampledFunction[m_?NumberQ, l_, n_?NumberQ, f_] :=
    If [ IntegerQ[Rationalize[n / m]], f, l ]

UpsampledFunction[m_List, l_, n_List, f_] :=
    Block [    {minv, mmatrix, nvector, retval},
        retval = l;
        nvector = Rationalize[n];
        If [ integerVectorQ[ nvector ],
             mmatrix = Rationalize[m];
             minv = If [ SameQ[mmatrix, lastMatrix],
                 lastMatrixInverse,
                 lastMatrix = mmatrix;
                 lastMatrixInverse = Inverse[mmatrix] ];
             If [ integerVectorQ[minv . nvector], retval = f ] ];
        retval ] /;
    Apply[ And, Map[ NumberQ, n ] ] && MatrixQ[m]

(*  UpsampleFactor --  must be independent of variable z  *)
UpsampleFactor[f_, z_, extractor_:GetAllExponents] :=
    Block [    {ratgcd},
        If [ FreeQ[f, z], Return[0] ];
        ratgcd = Apply[ RationalGCD, extractor[f, z] ];
        ratgcd = discBreakup[ratgcd, z];
        If [ IntegerQ[ratgcd],
             ratgcd,
             Numerator[ratgcd] ] ]

(*  UpsampleSequence  *)
UpsampleSequence[seq_, m_:2, fill_:0] :=
    Block [    {fillcount, input, len, newseq},

        len = Length[seq];
        newseq = List[ seq[[1]] ];

        (*  Build the new sequence newseq by appending
            m - 1 fill elements and then by appending
            a value from the original sequence.  *)

        For [ input = 2, input <= len, input++,
              For [ fillcount = 1, fillcount < m, fillcount++,
                AppendTo[newseq, fill] ];
              AppendTo[newseq, seq[[input]] ] ];

        newseq ]  /;
    IntegerQ[m] && ( m >= 2 )

(*  Z  *)
Z/: Z[n_] := Z[n, DummyVariables[Length[n], Global`z]]
Z/: Z[n_, z_] [ InvZ[zi_, n_] [f_] ] :=
    If [ SameQ[z, zi], f, f /. ReplaceWith[zi, z] ] /;
    Length[z] == Length[zi]

Z/: Z^-1[a__] := InvZ[a]

Z/: Format[ Z, TeXForm ] := "{\\cal Z}"
Z/: Format[ Z[n_, z_], TeXForm ] := StringForm[ "{\\cal Z}_{``}", n ]

(*  E N D     P A C K A G E  *)

End[]
EndPackage[]


(*  A L I A S E S  *)

Unprotect[ { AliasedSinc, AliasSinc, ASinc } ]
Clear[ AliasedSinc, AliasSinc, ASinc ]
AliasedSinc = Dirichlet
AliasedSinc::usage = Dirichlet::usage
AliasSinc = Dirichlet
AliasSinc::usage = Dirichlet::usage
ASinc = Dirichlet
ASinc::usage = Dirichlet::usage
Protect[ { AliasedSinc, AliasSinc, ASinc } ]

If [ TrueQ[ $VersionNumber >= 2.0 ],
     On[ General::spell ];
     On[ General::spell1 ] ]


(*  H E L P     I N F O R M A T I O N  *)

Block [    {newfuns},
    newfuns =
      { DeltaIntegrate,      DeltaPlot,        DiscreteGraphics,
        Extent1D,          DummyVariables,    GetDeltaFunctions,
        IntervalsToFunction,  OperatorInOperExpr,    OperatorName,
        ParametersInOperExpr, ScalingFactor,    SequencePlot,
        SequenceToFunction,   SignalPlot,        SignalsInOperExpr,
        "Signum",          SPSimplify,        TheFunction,
        ToContinuous,      ToDiscrete,        UpsampledFunction,
        UpsampleFactor,      UpsampleSequence };
    Combine[ SPfunctions, newfuns ];
    Apply[ Protect, newfuns ] ]

Block [    {newoperators},
    newoperators =
      { CConvolve,        CircularShift,
        Convolve,        DFT,              Difference,
        Downsample,        DTFT,            FT,
        Interleave,        InvDFT,            InvDTFT,
        InvFT,        InvL,            InvZ,
        L,            PolyphaseDownsample,    PolyphaseUpsample,
        Periodic,        Rev,              RewriteSummations,
        ScaleAxis,        Shift,              Summation,
        Upsample,        Z };
    Combine[ SPoperators, newoperators ];
    Apply[ Protect, newoperators ] ]

Block [    {newsignals},
    newsignals =
      { "AliasSinc",    "ASinc",        CPulse,
        CStep,        Delta,            Dirichlet,
        Impulse,        LineImpulse,        Pulse,
        Sinc,        Step,            Unit };
    Combine[ SPsignals, newsignals ];
    Apply[ Protect, newsignals ] ]



(* Under Mma 2.0, Delta is used by Integrate *)
If [ TrueQ[ $VersionNumber >= 2.0 ], Unprotect[Delta] ]

(*    When encoded using ":=", these next rules causes Simplify to   *)
(*  carry out incorrect simplifications when an expression contains  *)
(*  an Impulse term and a term containing 1/n because Simplify will  *)
(*  group the expression over a common denominator thereby putting   *)
(*  the n in front of Impulse[n] which whould force it to zero,      *)
(*  thereby eliminating a perfectly valid term.                 *)

Unprotect[ SPSimplificationRules ]
AppendTo[ SPSimplificationRules,
      ( f_ Impulse[n_Symbol + m_.] :> (f /. n -> -m) Impulse[n + m] /;
        NumberQ[m] && ! FreeQ[f, n] ) ]
AppendTo[ SPSimplificationRules, ( x_ Sinc[b_. x_] :> Sin[b x] / b ) ]
AppendTo[ SPSimplificationRules, ( Sinc[-x_] :> Sinc[x] ) ]
AppendTo[ SPSimplificationRules, ( Delta[-x_] :> Delta[x] ) ]
AppendTo[ SPSimplificationRules, ( Delta[a_?Positive x_] :> Delta[x] / a ) ]
AppendTo[ SPSimplificationRules, ( Impulse[-x_] :> Impulse[x] ) ]
AppendTo[ SPSimplificationRules, ( CPulse[L_, L_ - x_] :> CPulse[L, x] ) ]

(*  Simplifications when an impulse is the input      *)

AppendTo[ SPSimplificationRules,
      ( Downsample[L_, n_][a_. Impulse[n_]] :> a Impulse[n] ) ]
AppendTo[ SPSimplificationRules,
      ( Upsample[L_, n_][a_. Impulse[n_]] :> a Impulse[n] ) ]

(*  Simplifications when an exponential is the input  *)

AppendTo[ SPSimplificationRules,
      ( Rev[n_Symbol][Exp[ Complex[0, a_] b_. n_ ] ] :>
        Exp[-I a b n] ) ]
AppendTo[ SPSimplificationRules,
      ( ScaleAxis[l_,w_Symbol][Exp[Complex[0, a_] b_. w_]] :>
        Exp[I a b w l] ) ]

(*  Other simplifications:  when these are encoded    *)
(*  using ":=", they can cause rules in some of the   *)
(*  packages to confuse Mathematica to no end.  For   *)
(*  example, hard wiring these rules will cause       *)
(*  the "RewriteRules.m" package to hang.             *)
AppendTo[ SPSimplificationRules,
      ( Difference[k_, n_][c_] :> 0 /; MyFreeQ[c,n] ) ]
AppendTo[ SPSimplificationRules,
      ( Downsample[m_, n_][c_] :> c /; MyFreeQ[c, n] ) ]
AppendTo[ SPSimplificationRules,
      ( Periodic[k_,n_][x_] :> x /; MyFreeQ[x, n] ) ]
AppendTo[ SPSimplificationRules,
      ( Rev[n_][x_] :> x /; MyFreeQ[x, n] ) ]
AppendTo[ SPSimplificationRules,
      ( ScaleAxis[l_, w_][c_] :> c /; MyFreeQ[c, w] ) ]
AppendTo[ SPSimplificationRules,
      ( Shift[m_, n_][x_] :> x /; MyFreeQ[x, n] ) ]
AppendTo[ SPSimplificationRules,
      ( Summation[k_, rest___][a_ x_] :> a Summation[k, rest][x] /;
        FreeQ[a, k] ) ]

Protect[ SPSimplificationRules ]


(*  E N D     M E S S A G E  *)

Print["The knowledge base of signal processing operators has been loaded."]
Null